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ABSTRACT 



A procedure is presented that allows an engineer to coordinate 
protective devices consisting of fuses, circuit breakers, and over- 
current relays based upon load and fault currents. The resulting 
digital-computer program makes many of the decisions requiring very 
little prior in-depth knowledge of protective devices and power systems 
for its utilization. 

The coordination procedure consists of six separate programs 
that must run in sequence, as the output from one program forms the 
input for the next program. This reduces in size of internal computer 
storage required. 

The theory is presented with an explanation of each algorithm 
followed by a basic flow chart. An example computer run is presented 
at the end. 
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1. INTRODUCTION 



A prime objective of electrical engineers working on the design, 
operation, or maintenance of electrical power systems is to provide 
reliable, efficient power to their customers. The distribution system 
has always been a weak link in the power chain. Often, the catas- 
trophic results of a fault are easily seen, but the cause may be 
extremely difficult to find. A distribution system, either in the 
utility or the user's facilities that is not adequately protected from 
electrical faults is both dangerous and costly. An unsatisfactorily 
protected system can cause considerable inconvenience and damage to 
equipment. Adequate protection is not a set definition. Each load 
system is different with its own set of boundary values. 

In a distribution system, the flexibility of design reduces 
near the generator or the load. Usually the generator or load protec- 
tion is designed as an integral part of the equipment by the manufacturer 
and any changes should only be undertaken with their counsel. Moving 
away from the generator or load, the priority shifts from load pro- 
tection to distribution protection. It must never be forgotten that 
the load, distribution, and generation systems must be protected in 
such a way that a device that operates satisfactorily in one system 
cannot have a detrimental effect on another system. In this thesis 
the distribution section is attached from an infinite-bus power 
source to the load-protection device. The programs provide the 
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proper setting of a chosen device in order to provide good protec- 
tive coordination. Proper coordination is achieved when a device 
closest to the fault clears the fault or overload condition before 
the rest of the system is disturbed. In order to achieve this end by 
manual means, extensive mathematical calculations and experience are 
required. As a direct result, proper coordination is often not achieved 
The use of digital computers to solve the coordination problem makes 
a more exacting solution. However, computers with large internal 
memory are required to solve these problems. Several approaches 
to computer solutions have been published [1-4]. 

Problem Statement 

The purpose of this thesis is to develop the procedure and an 
interactive computer program for a generalized case of device coordina- 
tion. This coordination should include fuses, relays, and molded case 
circuit breakers. The phases of research are as follows. 

1. Develop a computer procedure to compute the one-line program 
per unit values. 

2. Develop a computer procedure to perform a load-flow analysis 

3. Develop a computer procedure to compute line-to-line ground 
and three-phase faults. 

4. Develop and apply procedures for mathematical representation 
of line X/R data, fuse, molded-case circuit breaker, and 
relay-curve time-current characteristics. 



I 
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5. Develop a computer procedure to coordinate the devices 
from phase 4 with the data derived from phases 2 and 3. 

Although the coordination programs are written for the PDP10 
computer, it requires no more than 25K bytes of internal memory and 
a disk unit to solve any system of 30 buses and 100 elements. Almost 
any large system can be reduced to its smaller parts composed of 30 
buses or less. As each segment is coordinated, it can be represented 
by a load and its protective device. In this manner a 300-bus system 
could be divided into smaller segments and in 11 runs would be com- 
pletely coordinated. The step-by-step approach used here allows for 
insertion or extraction of data before proceeding to the next step. 
Figure 1.1 is a basic flow chart of the coordination program. Here- 
after, it will be referred to by its individual parts. 

When the coordination process has been completed, a load flow 
study and a fault analysis has been performed as well. The major 
assumptions that have been made are listed here. Other minor ones 
are included in the corresponding chapters. These assumptions have 
been made in an effort to keep the solution as general as possible. 

1. The incoming line is assumed to be an infinite bus. 

2. Assymetrical currents are neglected. Devices are coordi- 
nated on the basis of symmetrical fault currents. 

3. Only synchronous-motor contributions to the fault 
currents are recognized. If a large induction motor 
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Figure 1.1. Flow chart of the coordination procedure. 
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Figure 1.1. 
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is present, it may be added as a fictitious synchronous 
motor. 

4. Line-charging shunt admittance is ignored. 

5. Prefault current and voltages are neglected in fault- 
current computations. 

6. Different sizes of fuses have the same time-current curve 
shape . 

7. Overload conditions are 125% of full load. 

8. The system has a balanced three-phase load. 

In this thesis each phase of research in the problem statement 
has a chapter devoted to it. After describing the computational 
methods used in curve smoothing, the coordination problem begins with 
data input. Next the fault currents are solved followed by the 
actual coordination program and device response graphing. Finally, 
a sample distribution system coordination problem is described. 



) 

s 




8 



2. CURVE SMOOTHING AND DIGITIZING 



General 



In this chapter, the methods for using a data digitizing "numonics H 
device and the curve smoothing program of Appendix A are explained. 

The various device equations that were found that adequately represented 
the particular device’s time current curve are also described. 

A method was described in Wagner [1] for performing this by 
using the PDP10 computer. This least-square curve-fitting program [1] 
was expanded to eliminate any external subroutines that may not be 
available. The coordination program makes extensive use of polynomial 
functions whose coefficients were derived by using the least-squares 
curve-smoothing technique . 

Curve Smoothing 

Electrical Cables . Coefficients for the mathematical modeling 
of resistance and reactance were taken from data for mine power 
cable [5]. This is used as only one type of cable or wire. Any other 
type may be modeled by the same procedure. For instance, either actual 
cable or line data could be used directly [6], or the theoretical resis- 
tance and reactance for each size could be calculated [5-8]. 
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Fuses . Two types of fuses are used in the coordination pro- 



gram. These are current-limiting and solid-material boric-acid power 
fuses. In each case, the curve shape is assumed to be the same for 
different sizes. The low-current asymptotic value is called the 
"of set." Each fuse has a minimum and a typical operating time. In 
fuse coordination, both values must be used in describing the operating 
envelope of a particular fuse. The higher typical operating time 
curve is referred to as M Ofset 1." By this means, any fuse can be 
completely described by four equations, where: 



Fuses are selected also by using mathematical models. To 
select a fuse, the program uses the overload-current value as the 
variable in the polynomial function. The result is the device number. 
In curve smoothing, variables can assume no value less or greater than 
those originally used as data for the curve-smoothing program, or large 
errors can result. 

Molded-Case Circuit Breakers. Like the fuses discussed above, 




n 

= ( E a. (Log Current -Ofset) J - 2 

j=° J 



( 2 . 1 ) 



n 



Log T - ( E a. (Log Current. -Of setl) 2 
yP k=0 



( 2 . 2 ) 



each molded-case circuit breaker consists of an operating envelope 
bounded by a minimum value and a maximum value. Unlike fuses, however, 
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a table look-up approach is used in selection of a device and its 
setting because of the smaller number of sizes and magnetic-trip settings. 

Relays . Wagner [1] discusses the representation of the family 
of relay curves for each type with one polynomial equation comprised 
of four coefficients. The 0^ power coefficient was left out as a 
variable for coordination. Unfortunately, this method proves unsatis- 
factory when working with more than one relay type. It was found that 

th 

each relay could be represented by five coefficients. The 0 power 
coefficient is necessary in that it provides the time displacement 
necessary for minimum time dial setting. The curves were derived 
using the procedure below. 

Westinghouse Company type relays were used. An average curve 
was taken as time-dial 4; This curve was reduced to a polynomial 
function. The time difference between time-dial 4 and time-dial 1/2 
taken at 20 times the pickup value is subtracted from the 0^ coef- 
ficient. Twenty times pickup was a good value for computing the time 
difference as the time value changes very little for increases in 
current beyond that point. The resultant curve was the average shape 
of a time-dial setting of 4 now shifted to minimum setting. This allows 
for time values computed during coordination to be summed with the 0^ 
coefficient. The curves are then represented as follows: 



Log Time = 



( a . (log(Current Value)^) + 
j=0 J 



Log (Time Setting) 
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This allows the use of more than one type of relay in each 
coordination as it separates time setting value from relay type. 

Curve Digitizing . The device’s time versus current curves 
were digitized by using a NUMONICS electronic graphics calculator 
(EGC) . The NUMONICS EGC is designed to translate graphic information 
on any medium into digital information. Complete information on its 
usage may be found in -NUMONICS Electronic Graphics Calculator [9] . 

The NUMONICS EGC then provided a convenient way to enter X-Y coordinate 
data representing any curve onto disk. Since this device works in 
length in inches, the data has to be normalized to its proper units. 

A simple program DIGIT. F4 was written to normalize and store 
onto disk the NUMONICS EGC information. The program flow chart is 
shown in Figure 2.1. 

Curve Smoothing . The curve-smoothing program contained in 
Wagner [1] was used as a basis for the curve smoothing program DEV.F4. 
Unnecessary parts were removed and a double precision matrix inver- 
sion routine was added. This prevented the necessity of using 
subroutines contained in the computer’s SSP routine files. The 
inversion routine was converted from The Pennsylvania State University 
Computer Center’s matrix inversion routine CMINV which inverts a 
complex square matrix. This routine in turn was converted from the 
SSP routine MINV [10]. 

Since it has its own inversion routine, this curve smoothing 
program can be easily expanded to any number of points and any 
dimension of fit by changing the dimension statements. The flow 
chart and procedure remains the same as in [1]. A complete program 
listing may be found in Appendix A. 
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Summary 



Polynomial functions describe the various device-response 
curves, line X/R data, and fuse-selection criteria. The coefficients 
for these polynomial functions are found from data entered using a 
digital-data device or manual input. This data is reduced by a curve- 
smoothing program to the desired function. 

The start of the coordination problem begins in the next chapter 
with line and load data entry. 
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3. THE INPUT OF DISTRIBUTION-SYSTEM DATA 



General 

Usually, the first step in performing any calculations on a 
distribution system is to reduce all system values to per unit. The 
methods used to do this in the program will be described in this 
chapter. Errors can be introduced easily when converting to the per 
unit system using manual calculations. Program INPUT. F4 was written 
to perform this function as well as computing line resistance and 
reactance values and storing all system data on disk for use in sub- 
sequent programs. 

Per-Unit Reduction 



The program uses the following formulas to convert to the per- 
unit system [8]. (See Table 3.1 for variable definitions used in this 
thesis . ) 



Z 



Z 



I 



Base 



P.U. 



Base 



(KV Base> /MVA Base 



z(n)/z 



Base 



KWA Base /CT Bas e 



) 
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Table 3. 


. 1. Variable definitions. 



z(n) 


element impedance 


z p.u. 


impedance in per unit 


x d 


sub transient reactance 


p 

p 


real power bus p 


% 


reactive power bus p 


E , 
P 


complex voltage bus p 


I 

P 


complex current bus p 


Y 

pq 


Y bus element bus p to q 


y 

pq,rs 


admittance value from primitive admittance matrix 


Z 

pq 


Z bus element bus p to q 


z pq>pq 


impedance value from primitive impedance matrix 


z f 


fault impedance 


y f 


fault admittance 
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KV . 2 KVA 

, Base given . . Base new . 

P . U . new Z P.U. given K KV^ ' rVA_ . ; 

Base new Base given 

Utilizing the above formulas it is possible then to represent all 
elements in the distribution system on a per-unit basis. Loads need 
to be included also in the one-line representation. 

The program allows loads to be represented either by induction 
motor horsepower, synchronous motor horsepower, kVA, or current. 

Loads are entered in the following format: Bus number/horse- 

power, kVA, current/power factor. The program assumes 100% efficiency 
so one must use: 



HP_ J = HP/Efficiency 

Entered J 



(3.1) 



If information on efficiency and power factors are unknown, the follow- 
ing equations may be used instead of 3.1 [9]. 

For induction motors, use: 



HP = HP x 1.1627 

Entered 



Power Factor = .86 



(3.2) 



For synchronous motors use: 

Unity Power Factor HP = 1.139 x Horsepower (3.5) 



.8 Power Factor HP = 1.18 x Horsepower 



(3.6) 
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Z 



P . U.new 



Z 



P . U . given 



KV . 2 KVA 

Base given . . Base new 

KV ' C KVA d 

Base new Base given 



Utilizing the above formulas it is possible then to represent all 
elements in the distribution system on a per-unit basis. Loads need 
to be included also in the one-line representation. 

The program allows loads to be represented either by induction 
motor horsepower, synchronous motor horsepower, kVA, or current. 

Loads are entered in the following format: Bus number/horse- 

power, kVA, current/power factor. The program assumes 100% efficiency 
so one must use: 



HP^ , = HP/Efficiency (3.1) 

Entered J 

If information on efficiency and power factors are unknown, the follow- 
ing equations may be used instead of 3.1 [9]. 

For induction motors, use: 



HP_ = HP x 1.1627 

Entered 



Power Factor = .86 



For synchronous motors use: 



Unity Power Factor HP = 1.139 x Horsepower 



(3.2) 



(3.5) 



.8 Power Factor HP = 1.18 x Horsepower 



(3.6) 
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Formulae for Line Resistance and Reactance 

Element data is entered in the following format: line bus/ 

load bus/cable length in feet/wire size/ transformer P.U. impedance/ 
size in M^A/line voltage in KV/load voltage in KV. 

Wire sizes that are entered in American wire gauge (AWG) are 
converted to MCM wire gauge by using the following formula [1] : 

2 

MCM = [(.005) [(1.229) (36 AWG) ] x 1000 (3.7) 



Once the wire size is in MCM, the line resistance and reactance are 
computed using the polynomial coefficients that were found by the 
procedure in Chapter 2. 



Resistance or reactance = Z a . [log (MCM) (3.8) 

j=0 J 



The program should not be limited to just one wire type. The program 
does this by reading wire coefficients a_, above from disk so that each 
time the program is run it can compute element data using different wire 
types until the entire system has been completely entered. 

Transformers become very important in that they can become the 
highest source of resistance and reactance. By their connections and 
neutral impedance, they determine the magnitude of zero-sequence current 
that will flow in the particular element. The program allows for this 
input and writes this information in the particular disk files for use 



later on in fault analysis. 
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Transformer inrush and withstand values must also be calcu- 
lated. These values must be used in coordinating protective devices. 
These devices must pass inrush currents but open under conditions of 
exceeded transformer withstand currents. The algorithm in the program 
is based on average values. If test data is available this should be 
entered directly by changing the transformer data disk file directly 
after the program is run and before proceeding with coordinations. 
Inrush and withstand values are computed as follows [11] : 



Inrush = 12x (Transformer Base Current) (3.9) 

Withstand @2 Sec = 25x (Transformer Base Current) (3.10) 

Withstand @5 Sec = 14 . 5x (Transformer Base Current) (3.11) 

Transformer withstand is some value between the two second value to the 
five second value and depends upon the transformer impedance. It is 
fairly easy, however, to have the device protect the complete range 
of withstand values and the program is written to accomplish this task. 
Synchronous subtransient reactance is computed as follows [11] : 

(Per unit KVA) x (.17) = iX, 

d 



(3.12) 



20 



One of the problems encountered in computing transmission line 
impedance is a reasonable calculation of the zero-sequence impedance. 
The computer program assumes that this value ± s two times 
the positive-sequence impedance [8] if nothing is entered when the 
computer asks for the zero-sequence impedance multiple. If a different 
value is entered, the computer multiplies the positive sequence reactive 
value by this multiple to arrive at the zero-sequence value. 

A basic flow chart follows in Figure 3.1 with a complete program 
listing in Appendix B. 

Summary 



The first step toward solving the coordination problem is com- 
pleted. All element data has been converted to the per-unit 
representation. The next logical step is to perform a load-flow 
analysis to determine the full load currents in each element. 
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Figure 3.1. 



Per-unit calculations flow chart. 
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Input 




Figure 3.1. Continued. 
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4. LOAD- FLOW ANALYSIS 



General 



In order to establish overload conditions on each element, 



the load currents must be determined. This is done by using a load 
flow analysis. The Newton- Raphson Method using Y-bus [12] was chosen 
because it required the least amount of iterations to obtain a solu- 
tion to the load-flow problem. 

Load Flow 



impedance network. When mutual coupling is ignored, Y-bus can be 
formed easily by algorithm. Any system can be completely represented 
by a set of equations [8,12]. 



In this method Y-bus must first be formed from the primitive- 




(4.1) 



E 



,bus 



and I 



bus 



are column matrices with one entry for each bus 



in the system. Y 



bus 



is a square matrix of order equal to the number 



of buses in the system. Therefore, once all bus voltages are known and 
Y-bus is formed, matrix multiplication gives the desired load currents. 
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However, there is no direct approach to obtain the bus voltages. 
The only approach to use is that of iteration. The Newton-Raphson 
method accelerates the iteration process so that normally an accurate 
answer is obtained after three iterations. 

The power at bus p is (from 4.1): 



P 

P 




* 

E I 
P P 



P 

P 




n 

E* I Y E 
P q=l Pq q 



(4.2) 



E 

P 



e +jf Y 

p p pq 



G — i B 

pq pq 



(A. 3) 



So 



P 

P 



n 

£ 

q = l 



{e (e G +f B ) + (f G -e B )} 
p q pq q pq q pq q pq 



(4.4) 



n 

~ ^ ( e G ~ff B )— e (f G -e B 'll 

P q = l P q pq q pq q pq q pq )} 



(A. 5) 



The Newton-Raphson method uses the Jacobian matrix of partial- 



differential equations as follows: 
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The elements of the Jacobian matrix are found from the following 
equations : 




e G -f B 
p pq p pq 



q^p 



9P n 

-r-2- = 2e G +f B -f B + I (e G +f B ) 

3e p p pp p pp p pp q pq q pq 

q=p 



(A. 7) 
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e B +f G 
p pq p pq 



q^p 



3P n 

ir^- = e B +Z G -e B + E (f G -e B ) 

3f p p pp fp pp p pp q=1 q pq q pq 

q^p 



(4.8) 
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Bus //I or the swing-bus voltage is known. The swing bus is 
the bus that will provide whatever power the system requires. The 
program begins by forming Y-bus by using the primitive-impedance matrix 
generated by INPUT. F4 of Chapter 3. Since mutual coupling is ignored, 
the y-admittance matrix can be found by inverting each member of the 
z-impedance matrix found in disk file TEST.DAT. y^ us matrix is formed 
using the algorithm below. 

Diagonal elements are formed as the sum of all admittances 
connected to that bus or node. Off-diagonal elements are the negative 
of the admittance connected between the buses or nodes. 

If more than 10 iterations are performed, the program stops and 
informs the user that either an error has been made in design or data 
entry. Upon successful completion of the program, a load-flow analysis 
may be obtained by printing file DISK.DAT. The program puts the load 
currents in the proper format and file for later use in the coordination 
program. A flow chart follows in Figure 4.1, and a complete listing 
is given in Appendix B. 

Summary 



This load-flow analysis obtained can be very useful in analyzing 
the distribution-system design. It can be used to point out locations 
for power factor correction as well as potential problems due to under- 
sized wiring. 

The load currents found here are stored on disk to be used in 
the coordination program. The next step before beginning the 
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coordination program is to solve for all fault currents. Z-bus must 
be formed to compute these currents. The methods used to do this 
follows in the next chapter. 




Figure 4.1. Newton-Raphson load-flow chart 
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3. FORMATION OF IMPEDANCE-BUS NETWORK 



General 



The formation of the impedance bus network (Z-Bus ) by 
algorithm is the first step toward solving for fault currents. In 
this chapter this algorithm will be discussed as well as other problems 
encountered in matrix compression. 

Formation of Z-Bus 



Once the one line diagram has been completed, the primitive 
impedance matrix that was found in INPUT. F4 is a partially filled 
matrix. However, it is not in the proper format for the formation of 
Z-bus as it contains no coupling. The primitive impedance matrix must 
change from one of order equal to the number of buses to one of order 
equal to the number of elements. This matrix becomes a very large 
sparse matrix. The memory required to store this matrix becomes too 
large to fit in a small computer. A method was derived to only store 
the non zero values and their location in the matrix. Since the matrix 
is triangular only the top half is stored (Figure 5.1). 

Diagonal elements found in INPUT. F4 are stored in "ZD." ZD 
is a three-dimensional matrix, ZD(200,2) where (X,l) is the element 
number and (X,2) is the diagonal value. 
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ZD(200 , 2) 




Number locations matrix (A) = (// elements) 2 X 2 

Number of locations of alternate method = 8 X (// elements) 



Figure 5.1. Transformation of 



square sparse primitive impedance. 
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The off diagonal mutual coupling Zm is a three dimensional 

matrix also where (X,l) is the location value A+jB. Element A is 

coupled to element B. Location (X,2) is the coupling value. Fortun- 
ately, partitioning makes it necessary to assemble a square matrix for 
inversion of only the coupled elements. Since coupling only exists in 
zero sequence, this applies only to the zero-sequence Z-bus matrix. 

In the formation of the Z-bus matrix, bus 1 is the reference 
bus for all calculations. Z-bus is formed by starting with element 
one and proceeding with each element in turn. An element is either a 

link or branch. A link means that both buses have been used earlier; 

a branch means they have not been used earlier. 

Like Y-bus in Chapter 4 any system can be represented by 



012 012 
Ebus = Zbus 



Ibus 



012 



Formation of Zbus can be formed either by the incidence, 
network-matrices method [12] or by algorithm. Formation by algorithm 
is much simpler and is used here. 

After forming the primitive impedance network, it is a straight- 
forward application of the equations of Table 5.1. 



[ 

1 
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Table 5.1. Equations for the formation of ZBUS . 

[ 11 ] 



Add p-q Mutual Coupling 



No Mutual Coupling 
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For example, a primitive impedance matrix looks like this 



for five elements: 



1(1-2) 2(2-3) 3(3-4) 4(1-4) 5(2-4) 



z 



o 



1 .l+j.5 j 5 

2 j 5 .l+j.5 

3 



.l+j.5 jlO j4 

j 10 -l+j.5 j8 

j 4 j 8 .l+j.5 



y 



O 




4 



5 



where the off diagonal values represent mutual coupling between ele- 
ments. Using this primitive matrix and the equations in Table 5.1 
yields Zbus. 



program asks for the zero-sequence impedance for all loads in the 
system and the mutual-coupling impedances. If the zero-sequence 
impedance of a load is not known, the program assigns it the same value 
as the positive sequence impedance. 



added are those of synchronous motors. The subtransient reactance 
is used instead of its actual impedance since the Z-bus is going to be 
used for fault calculations. Any link from a bus to the reference bus 
will cause a fault current. Prefault voltages and currents will be 
ignored as will all contributions from loads other than synchronous 
motors. A flow chart follows (Figure 5.2) and a complete program list- 



ZBUS.F4 is one of the longest programs in the series. The 



All loads that are added are links, and only loads that are 



ing in Appendix B. 




Zbus 



012 



Figure 5.2. 



formation flow chart. 
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Figure 5.2. 



Continued 
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{ Zbus^ 



formed} 



{All loads are links} 



Figure 5.2. Continued. 
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Summary 



Zbus is formed by algorithm allowing for mutual coupling 
between elements. The next step in the coordination procedure is to 
solve for fault currents. This procedure follows in the next chapter. 






] 

1 
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6 . FAULT CALCULATIONS 



General 



The final program required before executing the coordination 
program is FALT.F4. This program computes all the possible line-to- 
ground and three-phase- fault currents and prepares all other load 
data for entry into the coordination program. The procedure to 
accomplish this will be discussed in this chapter. 

Fault Currents . The equations describing the fault currents 
need to be derived [12] . 

Three-Phase Fault. For a fault at bus P: 



- - V 12 ( ™ w 01 S OI V\ (0) o 12 



( 6 . 1 ) 



where U is unity matrix, E (0) is prefault voltage, and Y is the 

P * 



fault admittance matrix for a three-phase fault. 



Y 



012 

F 




0 0 
1 0 
0 1 



( 6 . 2 ) 



AO 



Substituting (6.2) into (6.1) 



03.2 

L P(F) y F 



012 

1 P(F) y F 



012 

P(F) 



= 



0 

1 

0 



-1 



0 

0 

1 



1+Z 



PP l 
1+Z y. 
PP f 



1+Z 






0 

/T 

o 



0 0 0 
0 1 0 
0 0 1 



1+Z 



PP 



i+z 1 y T . 

PP f 



1+ z y r 

pp f 



o 

/ 3 
0 



0 

1 



i+Zpp yp 



0 

1 



i+z y., 

PP F 



0 

/ 3 
0 



(6.3) 



Since Zpp 



1 



Zpp and 



1/y f 



/ 3 



P(F) 



z f +z pp 









0 



(6.4) 
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Using the same method the following equations are found [12]. 



J P(f) 



J 3 2 F E PQ 
z F +Z PP 



(6.5) 



and 



i (F) 



= / 3 (E. 



i (0) 



z ip (1>E P «» 



z + z 

f pp 



( 1 ) 



i^P 




J p(f) 



°’ 2 = 0 



( 6 . 6 ) 



These equations are identical to those found for a three-phase 
grounded fault [12]. 

012 

The line-to-ground fault admittance matrix is given by Y 

r 

012 

below. Substituting Y into (6.1) the equations below result [12]. 

F 



012 = 

F 3 



1 1 
1 1 
1 1 



1 

1 

1 



012 



/ 3 E 



p(0) 



.0 



+2Z^+3z, 



PP PP 



(6.7) 



I 
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012 



/ 3 E 



Z < ‘°' ) +2Z ^ 1 ' > +3z 

PP PP F 



-z 


(0) 


PP 




+z 


(1) 


PP 




-z 


(1) 


PP 





( 6 . 8 ) 



r. 012 _ ^ 

E i(F) E i(0) 



i^P 



0 

0 



/ 3 E 



P(0) 



Z (0) +2Z (1) +3z 

PP PP -F 



Z. 

iP 

Z. 

iP 

Z. 

iP 



( 0 ) 

(1) 



(1) 



(6.9) 



The value of z is entered from the keyboard when asked by the program. 

r 

All of the equations are available now to solve for all fault currents. 
All post-fault bus voltages are solved. Using these bus voltages 
each element fault current is determined. These currents are stored 
onto disk so that they may be used in both the coordination and graph- 
ing programs. 

Summary 



The programs thus far have solved for all load and fault cur- 
rents but no coordination has been discussed. Now the coordination 
problem is ready to be addressed. In the chapter that follows, the 
protective devices for the distribution system will be coordinated. 



Start 
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Read Zbus 







^Set all bus 
voltages 


_\ 
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\ to disk 
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element 

currents 



\L 
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comp to three 
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Figure 6.1. Flow chart for FALT.F4. 
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7. PROTECTIVE-DEVICE COORDINATION 



General 



In this chapter, the protective-device-coordination problem 
will be discussed. A computational method is put forward that allows 
a digital computer to solve the problem. 

Device Coordination 



One of the best methods to approach coordination is a step by 
step device selection and setting working from the load toward the 
source. A load is the best place to start because there is not much 
flexibility for device setting at a load. The load flow analysis of 
Chapter 4 provides the overload values for the various loads. This 
enables a definite overload setting for each device protecting a load. 
Since overload currents are a primary concern at these buses, all fault 
settings can be placed at six times the full load currents. This should 
allow for normal starting currents. Care must be taken so that an 
instantaneous tripping element will not open under normal starting 
conditions. From the loads it is a straightforward process to place 
a line side device curve next to its load side curve to insure that 
proper coordination is achieved. By taking this approach, each device 
depends only on one other device. 



I 

1 
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This program is written with the option of entering two devices 
in each element — one located at each bus. From this approach, 
obviously two possibilities can occur. Either the device is on the 
line side of the element or on the load side of the element (Figure 7.1). 
If it is the line-side device, its setting depends only on the one 
load-side device. It must be coordinated for the fault current at 
the load bus and not the line bus [1] . Devices must be coordinated 
for a minimum possible fault current which will be discussed later. 

If the device is on the load side and it does not protect a 
load, then it must protect a following bus. This device must then be 
coordinated with all line-side devices connected to the bus it pro- 
tects. This load-side device must be coordinated for fault currents 
at the protected bus. 

Devices should be set to (1) open under overload conditions, 

(2) open under minimum fault conditions with proper coordination, and (3) 
open under maximum fault conditions with proper coordination. Often, 
however, all three conditions cannot be met simultaneously unless the 
proper devices are chosen. Criterion (1) is not as important as 
Criterion (2) or (3) when in the distribution system not at a load bus 
since an overload condition below fault conditions should not do damage 
to the wiring. For this reason wire sizes must be checked to insure 
that upon final coordination they can withstand the overload and fault 
conditions that can occur [5,13]. At this point it may be necessary 
to begin the entire coordination procedure again if an element cannot 
withstand the short-circuit conditions. Allowable short-circuit cur- 
rents can be checked by using the final coordination graph and the fault 



currents obtained in FALT.F4. 



I 

1 
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(a) Coordinating line side device of an element. 







c 

n i 


A 

O 


B 

r\ 


W 
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z 




0 1 
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1 






U 1 1 

Y 



(b) Coordinating load side device of an element feeding a bus with 
load elements connected W, Z. Y. 

Figure 7.1. Two possibilities for coordination element Z. In 

Figure (7.1a) coordination device A according to curve 
device B up to fault bus 2. In Figure (7.1b) coordina- 
tion device B depending on load flow current element Z. 
Coordinated with device curves C, D, E, up to fault at 
Bus 2. 
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This coordination routine utilizes fault currents and load 
currents in achieving protective-device coordination. The load cur- 
rents determine the overload or pickup setting in each device. The 
coordination problem is then in two parts. 

Devices must be coordinated for overload conditions. This 
provides the most important coordination step for fuses in that they 
contain no variable fault-protection settings as do relays or molded- 
case circuit breakers. Also, another reason is that based upon only 
fault currents and load-side device settings, it is impossible to 
satisfactorily achieve coordination. 

This inability to achieve coordination enters when setting a 
device that is to protect a bus with other elements attached to that 
bus as loads. Just laying a line curve next to a load curve with the 

highest time value at any current value will provide an erroneous 

answer. This device must be set by determining its overload current 
value first since this current is the sum of all currents in the load 
elements that are attached to that bus. Once this point has been 
established, each load device curve must be polled in turn to insure 
that the device setting responds to a fault in a load element slower 
than the load element line side device. The device needs to be 
coordinated up to the maximum fault current that the device will see. 

It is possible for this device to see a smaller fault current than that 
in a load element. 

Coordinating a device located at a line-side bus is much simpler 

since the pickup or overload setting is the same as the load-side 

device. The line-side device must be coordinated up to the maximum 
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fault current that will occur at the bus where its load-side device 
is located, A higher fault current will occur the closer one gets to 
the line-side device; however, this higher current value will only cause 
the line-side device protecting the element to respond quicker. The 
method of coordination will be to first find a load point based on a 
load flow, and second coordinate the device at any current value up 
to some maximum fault current such that the line-side device has a 
longer response time than its load-side device. 

Since the program coordinates three different types of devices, 
nine different possibilities can occur. Each device may be coordinated 
with one of its own type or one of the other two. 

The program structure lends itself well to the use of subrou- 
tines. Two of these are SELECT and SETDEV. SELECT chooses a device 
and sets its own load point, and SETDEV coordinates it. 

Fuse Selection and Coordination 



Due to the large number of fuse sizes available, look-up tables 
are not as efficient as using polynomial functions to select devices. 
Once a load current is known, this value is substituted into one func- 
tion providing the proper device number. This device number is next 
substituted into two other functions providing the OFSET and OFSET 1 
values referred to in Chapter 2. Since there are two types of fuses 
in the program, current-limiting and solid-material boric-acid power 
fuses, there are six functions used in selecting fuses. These all have 



the following form: 
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Y = I a.X j (7.1) 

3=0 J 

X will be the logarithm of overload current or device number. Table 7.1 

lists the various coefficients a.. 

J 

After the fuse is selected, it must be coordinated. The program 
uses the values of QFSET and QFSET 1 with the applicable fuse curve to 
solve for a time value for a set-current value. The program searches 
for the load-side device to find the highest time value to that current 
value. If the time value is larger for the size, it is properly coor- 
dinated. If not, it increases the fuse size by one, solves for QFSET 
and QFSET 1 , and tries again continuing until it has found the proper 
fuse size. 

Relay Selection and Coordinatio n 

Relays provide the most difficult coordination problem. First 
a turns ratio for the current transformer (CT) must be chosen. The 
computer selects the proper tap setting to prevent saturation. It 
does this by always selecting the highest turns possible while 
selecting the tap setting on the relay as well [12]. The computation 
uses a 600/5 CT [13] with tap settings and CO relays [14] with tap 
settings. This combination of tap setting for the CT and the relay 
combined with the load current establish the pickup point: 



Pickup =Log [(Turns Ratio CT x Turns Ratio Relay)] (7.2) 
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Table 7.1. Fuse-selection coefficients. 



Current to Device 



Offset Offset 1 



Current Limiting Fuses 






a o 


- 11.15221 


.9532 


1.1619 


a l 


23.865 


- .046433 


- .110148 


a 2 


- 16.487 


.117885 


.10626 


a 3 


2.9288 


- .029352 


- .02249 


a 4 


3.19672 


.00337 


.002452 




- 1.48877 


- .0001847 


- .0001365 


a 6 


.1776431 


.0000039 


.00000304 


Solid 


Material Boric Acid 






a o 


88.1643 


1.34359 


1.41734 


a l 


-215.2235 


.151393 


.1585889 


a 2 


188.225 


- .0239735 


- .035007 


a 3 


- 69.13403 


.00046627 


.007180 


a 4 


7.4695 


.0004321 


.0008126 


a 5 


1.65757 


.0000182 


- .0000377 


a 6 


.347786 


.0000078 


0 
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Since relay curves are in units of multiples of pickup this 
pickup value will be used often. Only the relay pickup has been set, 
its time dial must be coordinated. Table 7.2 lists the relay-curve 
coefficients. 

Due to the inverse time shape of a relay curve, many points on 
the curve must be checked with its load-side device to insure that 
proper coordination exists. The manner in which the family of relay 
curves are represented simplifies the problem somewhat. A relay 
operating time of 0.4 seconds is assumed throughout. At each current 
value, the load-side-device operating time is summed with the relay 
operating time. If this time is greater than the line-device time, the 
difference is the relay time-delay setting. Thus, when all currents 
up to the fault currents are checked, all points on the line-device 
curve are at least 0.4-seconds higher than all points on its load-side 
device . 

Molded Case Circuit Breaker 



The molded-case circuit breaker (MCCB) with its different 
characteristics provides a different set of problems. The MCCB curve 
cannot be represented by one equation. This device has one set of 
curves for the thermal element plus another for its magnetic element. 
The MCCB used in the program is a General Electric type, K-1200. 
Although there is some difference in the curve shape for the current 
ratings, 300-1200 A, they are similar enough to be represented as one 
type. Each MCCB has one curve corresponding to minimum total clearing 



time and a maximum total clearing time. 



Table 7.2. Relay-curve coefficients. 
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The magnetic portion of the time-current curve cannot be 
represented adequately by a polynomial function, but is represented 
by: 



If Current Magnetic Setting 

Minimum Time = 0, Maximum Time = .025 seconds. 

As will be seen later, .025-seconds response time is quite fast and 
usually does not present a problem during later coordination. 

The thermal elements of molded-case circuit breakers are 
susceptible to temperature change so the operating temperature must 
be entered. The temperature merely shifts the curve to the left or 
right so it does not affect the curve shape. This offset is added 
or subtracted from the multiple of current rating when using the thermal 
curve polynomial functions. 

The problem of coordinating the thermal elements is the same as 
that for fuses. The magnetic element coordination is much simpler in 
that the magnetic element is set separately. Its load-side device 
operating time is checked at a current value. If the operating time 
is greater, the magnetic element setting is increased. The end of the 
curve for the magnetic element is determined by the maximum fault cur- 
rent at the bus. The total decoupling of the overload from the fault 
device of a MCCB makes coordination much simpler than the fuse or relay. 
The more or less proportional- time shape of the curve will ensure that 
it will be properly coordinated once its proper current rating has been 
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established. Table 7.3 contains the polynomial coefficients used in 
the program. 

A basic flow chart (Figures 7.2 and 7.3) follows with a complete 
program listing in Appendix B. 

Summary 



In this chapter, the coordination problem was solved to allow 
fuses, molded-case circuit breakers, and relays to function so as to 
clear a fault in the area of the fault before interrupting the rest of 
the distribution system. There is enough flexibility in the program 
to allow device settings to be changed and checked using the program 
PK0TD.F4. An example problem would be beneficial at this point to 
show how to run the programs. In the next chapter, a sample distribu- 
tion system will be coordinated using the coordination programs. 
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Table 7.3. Molded-case circuit breaker. 



6 

log(time) = E a [log (current) - (ofset + . 1122) ] n - 2 
n=0 n 



6 

log(time) = E a . [log(current) - (ofset + .2315)^ -2 
3=0 J 



ofset = log(ckt breaker size) 



Minimum Time : 

a Q = 5.76908 

a x = -10.88261 

a 2 = 23.65418 

a 3 = -34.7922 

a, = 33.7747 

a 5 = -21.14536 

a, = 6.431778 

6 



Maximum Time: 



a Q = 5.8596721 

a = - 4.3943349 

a - - .50690283 

a 3 = -10.402739 

a, = -17.0649561 

a 5 = -13.011543 

a = - 4.0107239 
6 



COORD. F4 




device // 
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set device 




Input backup 
load protect. 
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Figure 7.2. Coordination flow chart. 
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Figure 7.2. 



Continued . 




Figure 7.3. Device-graphing flow chart. 
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8. OPERATING THE PROGRAMS 



General 



In this chapter the capabilities, limitations, and use of 
the programs will be described. The example system of Figure 8.1 
will be coordinated using these programs. 

One-Line Diagram 

INPUT. F4 . All data is entered via the keyboard in the format 
detailed in Chapter 3. Data entry is a straightforward process 
except when entering transformer data. The program is written for a 
"Y-Y M configuration. If one or both sides are "A" configured, then a 
large number must be entered for the transformer neutral impedance 
"(ZN)" on the bus connected to the A side. After the program is 
completed^ the zero sequence value for the transformer element con- 
taining the transformer may be assigned an arbitrary high value if the 
transformer wiring is such that there can be no ground current through 
the element. This is done by changing its value in disk file TEST.DAT 
(see Appendix C for the file format). 

When entering the loads the program asks which type of load 



to be entered. The various loads are described as: 
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Figure 8.1. One-line diagram of sample system. 
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HP - Synchronous motor 

HP1 - Assynchronous motor 
VA - In KVA 

CUR - In amperes 

S - No more loads 

Load Flow 



INPUT. F4 . There are no entries from keyboard for this program. 
Printing DSK.DAT will provide a load-flow analysis in per-unit values. 

Solution of Fault Currents 



ZBUS . F4 . In addition to forming Zbus , this program allows 
the input for mutual coupling and for zero-sequence impedance for loads. 
Mutual coupling tends to increase fault currents in the mutually coupled 
elements and small neutral-to-ground impedances can create higher line- 
to-ground fault currents than three-phase faults. 

FALT . F4 . Besides the solution of fault currents described in 
Chapter 6, FALT.F4 also allows the input of fault impedances. This can 
provide an interesting study of the system because the user can go 
back and place different fault impedances into the fault analysis to 
see if the protective devices will operate as desired. 



62 



Device Coordination and Plotting 



COORD . F4 . The program forms a two-dimensional complex matrix 
called ADEV that contains all the information required to plot all the 
devices in the system plus all actual settings that are to be made on 
the devices. An example row 3 in the matrix is 3, (3,4), (-4,2), 
(1.26,203), (1.03,1.26), (3,667,0), (.67,0), (20,0) in which: 

1. For 3, 3 is the row location of the matrix and corresponds 
to the element number. 

2. For (3,4), it is the first column location and is the bus 
numbers of the element. 

3. For (4, 2), -4 is a relay type // 4; 2 is a current limiting 
fuse type // 2. 

4. For (1.26,203), 1.26 is the logarithm of the pickup point; 
203 contains the tap information. Two is the tap of the 
CT, and 03 is the tap of the relay. 

5. For (1.03,1.26), 1.03 is the logarithm of the current for 
the minimum fuse open; 1.26 is the typical fuse open point. 

6. For (3.667,0), 3.667 is the magnetic setting of the relay; 

0 since the fuse had no magnetic setting. 

7. For (.67,0), .67 is the time-dial setting for the relay. 

The actual time dial setting must be read from the manu- 
facturer’s family of relay curves using the formula below 



at the maximum fault current. 
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time-dial-setting found = (found-time curve) - 
(time-dial 1/2 setting) 

8. For (20,0), 20 is the turns ratio of the current trans- 
former . 

Device numbers the programs asks for correspond to the 
devices below: 

1 - Current limiting fuse 

2 - Boric fuse high voltage only 

3 - Molded-case circuit breaker 

4 - Overcurrent relay 

In device 4 the program will ask later what particular type 
of relay to use from Table 8.1. 

Upon completion of coordination, a device listing is typed 
on the keyboard, and a printout of the same is available by typing: 
PRINT*. LPT. 

When performing the plot ; having a device listing as well as a 
fault current listing is helpful. 

PLOTD. F4 . This program plots the devices that were coordinated 
above. By using a graph, it is possible to visually check coordination. 
The program arrangement allows the user to change any device setting 
or size in ADEV.DAT to tailor the coordination to his particular require- 
ments (see Appendix C for the file format) . The program uses the lowest 
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Table 8.1. Overcurrent-relay types. 



Program Type // 


Relay Type 


i 


CO- 2 


2 


CO-3 


3 


CO-4 


4 


CO- 5 


5 


CO- 6 


6 


CO- 7 


7 


CO-11 



Relay type numbers are listed as negative numbers by the program. 
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voltage in the system as a base for plotting currents. The user can 
select any plot offset. This has the effect of plotting only the 
area of concern. In this way a very detailed graph can be obtained. 

Tables 8. 1-8. 4 provide a listing of program device type numbers 
to the actual device type and size. 

Table 8.5 details what is contained in the various disk loca- 
tions created by the programs. 

To coordinate the devices in Figure 8.1, the procedure begins 
by entering element data using INPUT. F4. The program asks if data is 
to be read from disk or from the keyboard. The first element must be 
entered using the keyboard. 

The proper procedure is to group the elements to be entered 
by wire or cable type as the program only handles one type of wire or 
cable at a time. When completed with one type, answer zero for all 
other questions and the computer writes the element data that it has 
computed on disk in the file RELA.DAT and stops. The new wire type 
polynomial coefficients are copied into RES. DAT and the program is 
ready to be run for the new wire type. 

The cable coefficients used in the test problem are for mine 
power cable [5]. The coefficients are listed in Table 8.6. 

After selecting keyboard, the base power and volts are entered. 

The computer then informs the user to begin element data entry. The 
data format (Chapter 3) is listed to aid the user. For example, the 
element from bus one to two of Figure 8.1 is entered as follows: 
1/2/3000/6. For the element containing the 1 MV A transformer the follow- 
ing is input: 2/3/35/6/ . 001 . 005/1/15/4 . 1. The bus base voltages are 



66 



Table 8.2. Current-limiting fuses. 



Program Type // 


Limiting Fuses 


i 


5E 


2 


7E 


3 


10E 


4 


15E 


5 


20E 


6 


25E 


7 


30E 


8 


40E 


9 


50E 


10 


65E 


11 


80E 


12 


100E 


13 


125E 


14 


150E 


15 


200E 
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Table 8.3. Boric fuses. 



Program Type // 


Size 


16 


15E 


17 


20E 


18 


25E 


19 


30E 


20 


40E 


21 


50E 


22 


65E 


23 


80E 


24 


100E 


25 


125E 


26 


150E 


27 


175E 


28 


200E 


29 


250E 


30 


300E 


31 


400E 


32 


2-250E 


33 


2-300E 


34 


2-400E 
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Table 8.4. Molded-case circuit breakers. 



Program // 



Size 



34 


300 


amps 


35 


350 


amps 


36 


400 


amps 


37 


450 


amps 


38 


500 


amps 


39 


600 


amps 


40 


700 


amps 


41 


800 


amps 


42 


1000 


amps 


43 


1200 


amps 
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Table 3.5. 



DIGIT. F4: 
DEV.F4: 
INPUT. F4 
LDFL0.F4: 
ZBUS.F4: 
FALT . F4 : 
COORD. F4: 
PLOTD.F4: 

RES. DAT: 
TEST . DAT : 
RELA.DAT: 
XFOR.DAT: 
LOCUR.DAT : 
DEV. DAT: 
DSK.XFR: 
DSK.DAT: 
ZBUS . DAT : 
ZBUS . FOR: 
MUT . CO : 
FALT . DAT : 
EFLT.DAT: 
LCUR.DAT : 
FCUR.DAT : 
ADEV.DAT: 



Disc files. 



Curve-smoothing program 
Numonics input program 
Element data-input program 
Load-flow program 

012 

Formation of ZBUS program 
Solution of fault-current program 
Coordinate device program 
Plot device curves 



Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 

Contains 



X/R line coefficients 

element and load data in per unit 

element and bus voltage data 

transformer inrush and withstand values 

load currents 

relay coefficients 

load-flow currents in P.U. 

load-flow analysis data 



ZBUS 



012 

012 



ZBUS~^~ before loads are added 
mutual coupling data 
per-unit fault currents 
per-unit fault voltages 
real values of load-flow currents 
real values of fault currents 
coordinated-device listing 
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Table 8.6. Mine power cable polynomial coefficients 



R-resistance 


X-reactance 



a Q 0 . 758690526 7 34D+01 


0 . 6 7314 52381407 621D-01 


a - . 1108155100578006D+02 


- . 1408237729509665D-01 


a 2 0.6290760668063740D+01 


- . 1721580927352306D-02 


a 3 - . 1627214281565102D+01 


. 6845126542245339D-03 


a. 0.1604746751660286D+00 

4 


0.0 
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input next with the voltage in volts. Upon completion of entry of all 
bus base voltages, the computer types the computed element value in 
per unit. The load data is entered next. The program types the data 
entry format for the user. The 50 hp motor is entered as: 8/50/.85. 

The only other data to be entered are mutual coupling and zero sequence 
load impedance in ZBUS.F4 and fault impedance in FALT.F4. The sample 
problem uses no coupling and a fault impedance of zero. COORD. F4 is 
run next. The device type desired is entered when asked for by the 
program. A sample data input is given in Appendix D. 
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The following is a list of load locations and 
currents . 



Bus number Load current 



4 84.0396680 
4 84.0396680 
7 293.6613400 
3 48.9435370 



The following is the current flow in the system 



under full 


load c ond i t 


ions as found from 


the 


load flow 


analysis. 








Bus P 


to Bus 


Q 


Current at P Current at Q 


1 




2 


34 .2226440 


34 .2226440 


2 




3 


34 .2226630 


125.2048700 


3 




4 


88.0372570 


88 .0372570 


3 




5 


39 . 17 227 10 


39 . 17 227 10 


5 




6 


39 . 172 1 2 1 0 


365.0129500 


6 




7 


313.4335700 


313.4835700 


6 




8 


51.6573760 


5 1 .6573760 



The following is the computed fault currents for a 
3 phase fault and a line to ground fault. 



Fit Bs P Element 3 Phase 



Ln to Gnd 



2 

2 

2 

2 

2 

2 

2 

3 

3 

3 



1 

2 

3 

4 

5 

6 
7 
1 
2 
3 



9380.7892000 
261 .6855700 
677.8190200 
280 .04 1 3400 
280.3767000 
2226 .3547000 
383.7966800 
1247.2941000 
1 247 . 2997000 
846.6528100 



7053.9298000 
183.7685400 
595.4521200 
187 .2686 100 
186.2585700 
1496.4606000 
270.9381000 
833 . 3620700 
326.2558200 
699.4013200 



4 

4 

4 

4 

4 

4 

5 

5 

5 

5 

5 

5 

5 

6 

6 

6 

6 

6 

6 

6 

7 

7 

7 

7 

7 

7 

7 

8 

8 

8 

8 

8 

8 

8 



-P> UJ K5 
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Klement 

4 

5 

6 
7 
1 
2 

3 

4 

5 

6 
7 
1 
2 

3 

4 

5 

6 
7 
1 



5 

6 
7 
1 
2 

3 

4 

5 

6 
7 
1 
2 

3 

4 

5 

6 
7 



3 Phase 
349.8265800 
349.7991 000 
2780.81 59000 
479.3659400 
1017.3596000 
1017.3546000 
4003.0938000 
285.3288800 
285 .3382800 
2268 . 1567000 
391 .0159800 
881 .8973500 
88 1 .8899900 
598.6094700 
3808 . 2506000 
353.1902700 
2807 .5742000 
483.9912400 
791.2624200 
791 .2645000 
537.0803000 
3416 .9009000 
3416.8793000 
2839 . 6040000 
489.5244300 
350.4617800 
350.4624900 
237 .3819000 
1 5 13 .3939000 
1513.3895000 
14273.8630000 
216.8161500 
156. 6019600 
156.6021300 
106.3009100 
676.2499600 
676 . 2485300 
561 .9951900 
6769 . 2755000 



Ln to Gnd 
206.0435600 
204 .2971400 
1648 .7648000 
302.5682200 
469 .6438400 
466 . 3491300 
2210 . 2807000 
121.9024600 
121 .3072000 
972.7699300 
174.0103800 
308.3295700 
307 . 2671700 
260.1660700 
1384 . 1725000 
108.5970700 
901 . 2339600 
183.8482000 
148.9287400 
148.4773900 
125.4932500 
668.5896000 
569.4760500 
517 . 3090900 
124 . 551 1 200 
103.1131400 
103. 1008700 
70.4041400 
445.7008900 
443.9770900 
4235.3970000 
64 . 2673200 
23.5717640 
23.5401980 
16.4766590 
102.2285800 
100.1874100 
84 .0780670 
1 059.3405000 
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Using the device listing of ADEV.DAT and the plotted curves, 
it is possible to obtain all device settings. 

Device number one is a molded-case circuit breaker. Its size 
is 300 A and its magnetic trip is on minimum setting. 

Device number 10 is a current-limiting fuse number 10. 

Device number 12 is a solid-material high-voltage boric fuse 
number 24. 

The relay data follows: 



Device 


Type 








Time 




Number 


Number 


Ct 


Tap 




Dial 


Magnetic 


3 


i 


3 (120) 


7 




1/2 


— 


4 


i 


3 


8 




8 


6 


5 


3 


3 


14 




5 


6 


6 


4 


3 


3 




4 


6 


7 


5 


3 


4 




6 


6 


8 


5 


3 


5 




6 


6 


9 


6 


3 


6 




10 


5.2 


2 


1 


3 


13 




5 


5.64 


11 


1 


3 


4 




10 


6 


13 


7 


3 


7 




2 


6 


14 


2 


3 


8 




1/2 


6 


The magnetic trip of device number 9 had 


to be reduced 


after 


coordination to allow 


the device 


12 to set 


properly to 


provide pro- 


tection 


for the 1 MVA 


transformer. The 


graphs that follow show the 


devices 


plotted to the 


lowest voltage in 


the system, 440 V. All the 


graphs were computer generated. 


The first 


graph shows 


how a one-line 


computer diagram would 


be placed 


to aid 


in 


understanding the graphs. 
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Lgure 8.2. Example coordination plot 



IME IN 
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gure 8.2. Continued 
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gure 8.2. Continued 



IME IN SECONDS 
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Summary 



The sample distribution system of Figure 8.1 was coordinated 
using the coordination programs. Using the device-plotting program 
PL0TD.F4 this coordination was checked to see if the computer generated 
coordination was satisfactory to the user. 

During the research for this thesis, many observations were made 
and many features could be added to these programs in the future. This 
will be discussed in the chapter that follows. 
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9. CONCLUSIONS 



In this thesis the author set out to represent a one-line 
distribution system in a per-unit system and to solve for three- 
phase and line-to-ground faults and to perform a load -flow analysis. 

The various protective devices-response curves needed to be represented 
by polynomial functions. These requirements were fulfilled using the 
methods detailed in Chapters 2-7. 

The last task was to develop a computer procedure to coor- 
dinate a distribution system’s protective devices. Fuses, molded- 
case circuit breakers, and relays were coordinated and their coordi- 
nation graphs were plotted. There are many extensions for further 
work that could be done. These extensions will be discussed later. 

While working with the various programs, several observations 
have proved useful when performing coordination. Coordination should 
not be an afterthought in designing a distribution system. It is 
easy to design a system that is impossible to coordinate correctly. 

A radial system is easier to coordinate. Moving from a load toward 
the source, the time-current curves shift to higher current-time 
responses so that when finally arriving at a transformer, providing 
proper transformer protection and coordination simultaneously may be 
impossible. 

Coordination of the thermal-tripping elements of molded case 
circuit breakers is very difficult due to the difference between the 
minimum and maximum tripping time curves. Coordination of these 
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devices is best done manually with some overlap of the thermal time 
curves. The magnetic elements could then be adjusted to provide 
proper fault coordination. Those manually adjusted devices can be 
added to ADEV.DAT and be checked using PL0T.F4. 

In systems with different voltages, the transformers should be 
of different sizes so that the smaller sizes are on the load side. 
This will allow proper coordination and protection. Relays are 
preferable to other types of protection on the high-voltage parts 
of the system in that their response can be controlled. 

Recommendations for Further Work 



1. Using an analog computer and the programs presented here, 
an entire distribution and generation system could be 
looked at in detail. A transient-stability study would 
show the effect of various device settings. 

2. The INPUT. F4 program could be enlarged to automatically 
account for different transformer configurations. 

3. Additional devices could be added easily to the program. 

The program allows up to twenty different relay types. The 
additional coefficient data is merely entered in disk 

DEV. DAT. Additional molded-case circuit breakers would be 
slightly more difficult in that it would require some program 



modification . 
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4. The plotting program could be enlarged to include drawing 
a one-line diagram in the space provided and adding other 
data directly on the graph. 

5. The program could be expanded for the case of having one 
device in an element. 
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CURVE -SMOOTHING PROGRAMS 
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C 

* * * * * * * 

l 

5 

8 

10 

15 

20 

30 

35 



100 



150 

200 

210 
22 0 



TITLE : DIGIT. F4 ; CREATES FILES FROM DATA 
USING HUMONICS ARM 

********************************************* 
TYPE 5 

FORMAT ( ' INPUT DEVICE NR..') 

ACCEPT 8 , I D E V 
F0RMAT( I ) 

OPEN(UNIT=IDEV) 

TYPE 10 

FORMAT ( ' INPUT X-AXIS DIMENSION, Y-AXI S 
DIMENSION' //'..,') 

ACCEPT 1 5 , XAX , Y AX 

FORMAT ( 2 F ) 

TYPE 20 

FORMAT ( ' INPUT LENGTH PER LOG 1 0 ; X , Y . . ' ) 
ACCEPT 35 , XL , YL 
FORMAT( 2 F ) 

XAX=100*XAX/XL 

YAX=100*YAX/YL 

CALL 0PEN(30,' DEVICE',' NUM' ,' MODE ',' IMAGE ' , 
l'ACCESS'.'SEQ IN', 'BUFFERS', 4) 

CALL NUMRES( 30,4) 

CALL LOCK( 1,1) 

1< E A D ( 3 0 , E N D = 2 0 0 ) IDATA 

I X=( I DATA/" 100 000 0) .AND. " 37 777 

IY=IDATA. AND. "37777 

IF LAGS=( IDATA. AND. "7 4 0000)/ "40000 

X=IX/XAX 

Y=I Y/YAX 

WRITE ( I DEV, 150) I FLAGS , X, Y 
FORMAT ( I , 2 F ) 

GO TO 100 
CALL NUM CL 0 
CALL U N L 0 C K ( 1,1) 

CLOSE( UNIT=IDEV) 

TYPE 210 

FORMATC' CONTINUE? ..,'$ ) 

ACCEPT 220, IC 
FORMAT( I ) 

I F ( IC.EQ.O)GO TO 1 

STOP 

END 
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ft***********************************************}!?*** 

C CURVE FITTING PROGRAM : DEV . FA 

;fr********:k**7W< *&*&&************* ******************** 

INTEGER L( 7 ) , MM( 7 ) 

REAL X( 100) , YS( 100) 

REAL* 8 YM(100),G(7,7),U(7),A(0/6),GG(7,7), ERROR ( 100) 
TYPE 100 

100 FORMAT( 1 X, ' CHOOSE DIMENSION OF FIT '$) 

ACCEPT 1 1 0 , N D 
110 FORMAT ( I ) 

TYPE 120 

120 FORMAT( IX, ' CHOOSE NO. OF POINTS '$) 

ACCEPT 1 1 0 , M 
TYPE 130 

130 FORMAT( /' CHOOSE INPUT FILES'/' •••, '$) 

ACCEPT 1 1 0 , I UN I T 1 
DO 140 I = 1,M 

RE A D ( I U N I T 1 , 135)J,X(I) , Y S ( I ) 

135 FORMAT ( 1 , 2 F ) 

Y S ( I ) = ALOG 1 0 ( YS ( I ) ) 

140 CONTINUE 

DO 10 K= 1 , M 
DO 10 1 = 1 , N D 

U(I)=U(I)+YS(K)*X(K)**(I-1) 

DO 10 J= 1 , ND 

C( I , J ) = C( I , J ) + X ( K ) * * ( I+J-2 ) 

10 CONTINUE 

CALL DMINV ( G , 7 , ND) 

DO 20 1=0 , ND- 1 

DO 20 J = 1 , N D 

A( I )=A( I)+G(I + 1 , J ) * U ( J ) 

20 CONTINUE 

DO 30 1 = 1 , M 
DO 30 J = 0 , ND- 1 

Y M ( I) = YM( I ) + A(j )*X( I )**J 

30 CONTINUE 

DO 4 0 1 = 1 , M 

40 SQUERR=SQUERR+( YS( I )-YM( I ) )**2 

TYPE 160, SQUERR 

160 FORMAT( ' THE I NTERGRAL- SQUARED ERROR IS: ',E) 

OPEN(UNIT=40,FILE='A.DAT' ) 

WR I TE ( 4 0 , 170)A 
CLOSE(UNIT=40) 

170 FORMAT ( D ) 

TYPE 99 , IUNIT1 
99 FORMAT( ' IUNIT=' , I) 

DO 50 1 = 1 , M 

5 0 ERROR( I ) = DABS( (YM( I ) — YS( I) ) / Y S ( I) ) 

OPEN( UN1T=4 1 , FILE=' YINOUT. DAT' ) 

WRITE (4 1 , 180)(X(1) , Y S ( I ) , Y M ( I) , ERROR (I) , 1 = 1 , M ) 

180 FORM AT ( 2 F , 2D) 
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C 

MATRIX 

C 



C 

C 

C 



10 

15 

20 



25 



30 

35 

38 



40 

45 

50 

55 



60 
6 2 
65 



7 0 
75 



CL0SE( UNIT=4 1 ) 

STOP 

END 

TITLE : DMINV ; SUBROUTINE TO INVERT A REAL SQUARE 
DOUBLE 
PRECISION 

SUBROUTINE DM I NV ( A , MD I M , N ) 

RE AL*8 A( MDIM , N) 

INTEGER L( 100) , M( 100) 

CONVERTED FROM SSP ROUTNE MINV 

DOUBLE PRECISION BIGA , HOLD 
DO 80 K= 1 , N 
L ( K ) = K 
M ( K ) = K 
BIGA=A( K, K) 

DO 20 J=K,N 
DO 20 I = K , N 

IF(DABS( BIGA)-DABS(A( I , J) ) ) 15,20,20 

BIGA=A( I , J ) 

L(K)=I 
M( K)= J 
CONTINUE 
J=L(K) 

IF(J-K) 35,35,25 
DO 30 1 = 1, N 
HOLD=-A(K, I ) 

A(K,I)=A(J,I) 

A( J , I ) = HOLD 
I = M ( K ) 

IF(I-K) 45,45,38 
DO 40 J = 1 , N 
110LD = -A( J , K) 

A( J , K) = A( J , I ) 

A( J , I ) = H 0 L D 

DO 55 1=1, N 

IF(I-K) 50,55,50 

A( I , K)=A( I , K)/(-BIGA) 

CONTINUE 
DO 65 1=1 , N 
!10LD = A( I , K ) 

DO 65 J=1 , N 
IF(I-K) 60,65,60 
IF(J-K) 62,65,62 
A(I, J)=H0LD*A(K, J)+A(I , J) 

CONTINUE 
DO 7 5 J = 1 , N 
IF(J-K) 70,75,70 
A(K, J)=A(K, J)/BIGA 
CONTINUE 



8 0 


A(K,K)=(1.0000000)/BICA 

CONTINUE 


1 00 


K= N 
K= K- 1 


105 


I F ( K ) 1 50, 150, 105 
I=L(K) 


108 


IF(I-K) 120,120,108 
DO 1 10 J= 1 , N 


1 1 0 


HOLD=A( J , K) 

A( J,K)=-A( J, I) 
A( J , 1 )=HOLD 


120 


J = M ( K) 


125 


IF(J-K) 100,100,125 
DO 130 1 = 1, N 


1 30 


H0LD=A( K, 1 ) 
A(K, I)=-A( J, I) 
A( J, I ) = HOLD 


1 50 


GO TO 100 
RETURN 




END 



APPENDIX B 



COORDINATION AND PLOTTING PROGRAMS 
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***************** *******************&****** ************ 

C 1NPUT.IT inputs data from a terminal, places in 

C THE PROPER FORMAT FOR OTHER PROGRAMS 

******************************************************* 

c 

c 

DIMENSION 

L ( 3 0 , 30,2) , P U ( 3 0 , 3 0 , 1) ,RES(30,30,2),E(30) ,XL( 3 0,30,2) 
DIMENSION 

S(300) , ZER( 30, 30, 2) , RC(7 ) ,XC( 7 ) , ZG(30) ,XD(30) 

COMPLEX T1,TPU,PU,S,S1,PU1,PZ,CUK,ZER,ZG 

REALMS A ( 7 ) , AA( 7 ) 

REAL L 

********************************************************** 

C IF THE PROGRAM HAS ALREADY BEEN RUN ONCE ALL DATA 

C NEED NOT BE ENTERED AGAIN 

*********************************************************** 

TYPE 3 

3 FORMAT( ' DO YOU WANT TO READ DATA FROM A FILE OR 

FROM KYBRD? ' / 

2' < C R > F 0 R INPUT FROM TTY , 1 FOR READ...') 

ACCEPT 11,1 
IF(I.NE.0)GO TO A 00 

************************************************************ 

C INPUT OF DATA ON LINES 

******************************************************** 

TYPE 5 



FORMAT( 2 X , ' ENTER 
ACCEPT 10, ZBASE 
TYPE 6 


THE 


PROBLEM 


BASE 


IN 


MV A' ) 


FORM AT ( 2 X , ' ENTER 
ACCEPT 1 0 , E B A S E 


THE 


PROBLEM 


BASE 


IN 


KV' ) 



10 FORMAT(F) 

11 FORMAT(I) 

GO TO 14 

***************************************** * *************** 

C PROGRAM RETURNS HERE AFTER READING DATA AT 400 

********************************************* * * * * -A ****** 

12 TYPE 13 

13 FORMAT( ' OLD DATA IS ENTERED , READY FOR NEW LINE 
DATA. . ' ) 

ACCEPT 11,1 

14 TYPE 15 

15 FORMAT( /, 2X, ' ENTER ELEMENT DATA AS 

FOLLOWS : BUS( )TO( ) ( LENGTH) 

1 ( SP) S IZE' , / , 2X, ' (TRANSFORMER 
I M P E D E H C E ( P.U. ) ) (TRANSFORMER SIZE 

2 1 N MV A) SP',/, 2X, ' (TRANSFORMER HIGH VOLTAGE IN KV) 

SP' , / , ' 

3( TRANSFORMER LOW VOLTAGE IN KV ]...,'/ ) 

’k'k'k'k'k'k’k’k'k'k’k It 'k'k'k'k'k'k'k’k'k'k'k'k'k'k’k'k’k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k’k’k’k’k'k'k 
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C THESE TWO FILES CONTAIN COEFFICIENT DATA FO AO BELOW 

C THEY ARE CREATED USING THE CURVE FITTING 

PROGRAM : CURV . FA 

****** **************************************************** 
OPEN( UNIT=A0 , FILE=' RES. DAT' ) 

OPEN(UNIT=35 , FILE='XFOR. DAT' ) 

READ( AO, 20) ( A( I ) , 1 = 1 , 7 ) 

READ( AO, 20) ( AA( I) , 1=1 , 7 ) 

20 FORMAT(D) 

DO 21 1=1,7 

RC( I)-A(I) 

XC( I ) = AA ( I) 

21 CONTINUE 
CLOSE( UNIT=A0) 

C L 0 S E ( U N I T = A 1 ) 

22 ACCEPT 2A , I , J , XLL, WIR, 1 1 , TS] , TVH, TVL 

2 A FORMAT( 2 1 , 7 F) 

IF( I . EQ. 0) GO TO 1 10 

TYPE 26 

JJ=1 

26 FORM AT ( ' PARALLEL LINE 1 TO 3...,',$) 

ACCEPT 27, JJ 

27 FORMAT(I) 

IF( JJ.EQ.O) JJ=1 
L( I , J , J J ) =XLL 
IF( WIR .EQ. 1 AO) GO TO 31 
I F ( WIR. EQ. 1 30 ) G 0 TO 32 
IF( WIR.EQ. 120)GO TO 33 
IF( WIR.EQ. 1 10)CO TO 3A 
S A = W I R 
S B = W I R 

IF( SA.LT. 37) GO TO 35 
GO TO AO 

31 SA=-3 

GO TO 35 

32 S A= - 2 

GO TO 35 

33 SA=— 1 

GO TO 35 
3 A S A = 0 

GO TO 35 

****************** ********* ************ ************** 

C CONVERSION FROM AWC TO MCM 

***************************************************** 

35 SC=1 . 1229**(36-SA) 

SD= . 005*SC 
SE=SD**2 
S B= 1 000* SE 

******************************************************** 

C FORMULAE FOR LINE RESISTANCE AND REACTANCE 

******************************************************* 
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AO 

EES ( I , J , J J.) = RC( 1 ) + RC( 2 )*( ALOG10( S B) ) + RC( 3 ) * ( ALOC 1 0 ( S B 
)** 2 )+ 

1 RC( 4 )*( ALOGIOC SB) **3 ) + RC( 5 ) *( ALOG1 0( SB)**4 ) + RC( 6) *( A 
LOCI 0 (SB) 

2**5 )+RC(7 ) * (AL0G10( SB) **6 ) 

XL( I , J , JJ)=XC( 1 )+XC( 2 )*( ALOCIO(SB) )+XC( 3 ) * ( ALOG 1 0 ( S B ) 
** 2 ) + 



1 XC(4 )*( ALOG1 0( SB) **3 ) + XC( 5 )*( ALOGIOC SB) **4 ) + RC( 6 ) * ( A 
LOG 1 0 ( S B ) 

2 * * 5 ) + X C ( 7 ) * ( ALOG 1 0 ( S B ) * * 6 ) 

TYPE 50 , RE5( I , J, JJ) , SA, SB , XL( I , J , JJ ) 

50 FORMAK 2X, ' LINE RES. ' , 3F,/,'LINE REACT.. ' ,F) 

I F ( CABS(T1 ) . EQ . 0)GO TO 85 

PU(I,J,JJ)=T1*(ZBASE/TS1)*((TVH/EBASE)**2) 

TYPE 55, I, J 

55 FORMAT(' INPUT TRANSFORMER ZN IN P.U. REFER TO BUS 

SIDE, ( BUS 

1 ' , 13 , ' ) SP( BUS ' , 13 ) 

ACCEPT 6 0 , S 1 , P U 1 
60 FORMAT ( A F) 

TV = TV H* 1 E3 

ZG( I) = S1*(ZBASE/TS1 )*( ( EBASE/ TVH ) ** 2 ) 

ZC( J)=PU1*(ZBASE/TS1 )*( ( E B AS E/ T VL ) * * 2 ) 

CURB=TS 1 * 1 E6 /TV 
TINRSI1=1 2*CURB 
TWST1=25*CURB 
TW5T2= 1 A . A*CURB 

WRITE (35, 70) I , J, TINRSH, TWST1 , TWST2 , TS1 
70 FORMAT( 21 , AF) 

85 IF( IMAX.GT. I ) GO TO 90 

IHAX= I 

90 IF( JMAX.GT. J) CO TO 100 

J M A X = J 

100 ZER( I , J , JJ)=CMPLX(0 . 0 , XL( I , J , JJ) ) 

TYPE 102 ,ZER( I , J, JJ) , SSI 
102 FORMAT ( 3F) 

GO TO 22 

110 0 P E N ( U N I T = 3 1 , F I L E = ' RELA.DAT' ) 

OPEN(UNIT=32 , FILE='TEST. DAT' ) 

OP KH( UNIT=33 , FI LE= ' LOCUR . DAT ' ) 

IF( IMAX.LT. JMAX) GO TO 115 
I B = I M AX 
CO TO 118 
115 I B= JMAX 

C BUS VOLTAGES ARE INPUT 

*******************&******>£************* ***************** 
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113 DO 120 1=1, IB 

TYPE 119,1 

119 FORMAT ( INPUT BUS ', I 3 ,' VOLTAGE ...,',$ ) 

ACCEPT 10, E2 

I F ( E 2 . EQ . 0 ) GO TO 121 
E ( I ) = E 2 

120 CONTINUE 

121 TYPE 122 

122 FORMAT( / , 2X, ' DATA IS ACCEPTED AND A LIST FOLLOWS') 
TYPE 1 2 5 , ZBASE 

125 FORMAT( 2X, ' PROBLEM BASE IN MVA..,',F) 

WRITE (31, 126)ZBASE,EBASE, IB 

126 FORMAT( 2 F , I ) 

El = E( 1 )/(EBASE*lE3) 

WRITE ( 32 , 1 28) , IB , El 
128 FORMAT( I , 2F) 

DO 135 1=1,30 
DO 133 J= 1 ,30 
DO 133 J J= 1 , 2 

IF(L(I, J, JJ).EQ.0)CO TO 133 
NN=NN+1 

P Z = ( C MP L X ( R E S ( I ,J,JJ) , XL ( I , J , J J ) ) ) * ( .001*L(I,J,JJ)) 
PZZ=( ZBASE* 1E6)/(E( I)**2) 

PZ=PZ*PZZ 
T1=PZ+PU( I , J , JJ) 

ZER( I , J , JJ)=ZER( I , J , J J)*PZZ*( .001*L( I , J , JJ) )+Tl 
TYPE 130,NN,L(I, J, JJ) ,T1 

130 FORMAT( 2X, ' ELEMENT ' , I 3 , 5 X , ' LENGTH ' , F , / , 

1 IX, ' TRANSFORMERPER UNIT IMPEDENCE PLUS LINE 
IMPEDENCE' , 2F) 

WRITE (31, 131)1, J,JJ,L(I,J,JJ),RES(I,J,JJ),PU(I,J,JJ), 

XL(I, J, JJ) 

1 , ZG ( I ) , ZC ( J ) 

131 FORMAT( 3 I , 9F) 

WR I T E ( 3 2 , 1 3 2 ) NN , I , J , J J , T 1 , Z ER ( I , J , J J ) , Z G ( I ) , Z G ( J ) 

132 FORM AT ( 4 I , 8 F ) 

133 CONTINUE 

135 CONTINUE 

W R I T E ( 3 1,136) 

WRITE( 32 , 136) 

136 FORMAT( ' 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0') 

DO 138 1 = 1, IB 

WRI TE( 3 1 , 137)1, E(I) 

137 FORMAT(I,2F) 

138 CONTINUE 

WRI TE ( 3 1 ,139) 

139 FORMAT( ' 0 0.0 0.0') 

TYPE 1 4 0 , I B 

140 FORMAT( 1 OX ,' NUMBER OF BUSSES IS...,', I) 

145 TYPE 150 
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&****£*&**************£********************&&£****£****** 
150 FORMAT ( / / ' NOW ENTER THE LOAD DATA IN THE FOLLOWING 

MANNER; 

l'/' (BUS ) ( LOAD IN K V A ) ( P W R FACTOR) ... OR .. . 

2'/' (BUS ) (LOAD IN HP)(PWR FACTOR). OR....'/ 

3' (BUS ) ( LOAD IN AMPS)(PWR FACTOR)') 

DO 300 N = 1 , I B 

TYPE 160 

160 FORMAT ( ' WHAT METHOD WILL YOU USE 

? . .VA,HP,OR AMPS? ? ' ) 

ACCEPT 162, MD 
162 F0RMAT(A4) 

TYPE 164, MD 

164 FORMAT( 2X, A4 , ' ENTER...') 

I F ( MD . EQ . ' HP1 ' ) GO TO 170 
IF( MD. EQ. ' HP ' ) GO TO 170 
IF( MD. EQ. ' AMPS' ) GO TO 180 
IF( MD . EQ . ' S' ) GO TO 305 
MD= ' VA' 

ACCEPT 165, J , AS 1 , PF AC 

165 FORMAT ( I , 2F) 

SZTZ=AS1/(ZBASE*1E3) 

SZR=SZTZ*PFAC 
PCOS=ACOS( PFAC) 

SZI=SZTZ*SIN(PCOS) 

S( J )=CMPLX( SZR, SZI) 
AI=AS1*1E3/(E(J)*SQRT(3.0)) 

WRITE(33 , 168) J, AI 

168 FORMAT ( I , F ) 

ST=MD 
GO TO 300 

170 ACCEPT 175 , I , PRH, PFAC 

175 F0RMAT(I,3F) 

SZ TZ = PRH* 7 4 6 
SZ R= SZ TZ * PFAC 
PCOS=ACOS( PFAC) 

SZI=SZTZ*SIN( PCOS) 

S(I)=CMPLX(SZR, SZI)/ (ZBASE*1E6) 

IF( MD. EQ . ' HP' )XD( I)=. 17/CABS (S( I) ) 

AI=SZTZ/ ( SQRT( 3 .0)*E( I ) ) 

WRITE( 33 , 178) I , AI 
178 FORMAT ( I , F ) 

ST = MD 
GO TO 300 

180 ACCEPT 185 , I , CURT, PFAC 

185 FORMAT( 1 , 3 F ) 

WRITE( 3 3 , 187 ) I , CURT 
187 FORMAT(I,F) 

CURll=CURT* PF AC 
PCOS=ACOS( PFAC) 

CURI=CURT*SIN( PCOS ) 



300 

305 



3 10 

312 

315 

320 



330 

331 

332 
335 
400 



4 1 1 
4 50 
4 56 



470 

472 

457 

458 



9 00 



CUR=CMPLX( CURR, CUR I ) 
CUR=CUK*(SQRT(3.0)*EBASE/ZBASE) 

S( I )=3*EBASE*CONJG( CUR) / SQRT( 3 . 0 ) 
S T = M D 

CONTINUE 
DO 320 N = 1 , I B 
NN = N 

TYPE 310, N, S(N) , MD 
FORMAT( ' BUS ' , I 3 , 2 F , A4 , 2 F ) 

WRITE ( 3 1 , 31 2)NN, S(NH) 

FORMAT( I, 2F) 

WRITE( 32 , 3 1 5)N, S(N) , XD(NN) 

FORMAT( I , 3 F ) 

CONTINUE 
WRI TE (31 ,330) 

WRITE ( 32,331) 

FOR.MAT ( ' 0 0.0 0.0') 

FORMAT( ' 0 0.0 0.0') 

WR I TE ( 3 3 , 332 ) 

FORMAT( ' 0 0.0') 

W R I T E ( 3 5 , 335 ) 

FORMAT( 0 0 0.0 0.0 0.0 0.0 ') 

CO TO 900 

OPEN( UNIT = 3 6 , FILE= ' RE LA . DAT' ) 

READ( 36 , 4 1 1 ) Z BASE , EBASE , I B 
IMAX=I B 
F ORMAT ( 2 F , I ) 

READ( 3 6 , 456)1, J,JJ,XL3 ,RES1 ,PU1 , XL 2 , T 1 , S 1 
FORMAT (31 , 9F) 

IF(I.EQ.O)GO TO 470 

L( I , .1 , JJ ) = XL3 

RE S ( I , J , J J ) = RES1 

PU ( I , J , 1 )=PU1 

XL ( I , J , J J ) =X L2 

ZER(I,J,JJ)=CMPLX(0.0,XL2) 

ZC( I ) =T1 
ZG ( J ) = S 1 

IF( I . NE . 0 ) GO TO 450 
READ (36, 472)1, E(I) 

FORMAT( I , 2F ) 

I F ( I.NE.O) GO TO 470 
READ( 36 , 458 ) I , SI 
FORMAT( I , 2 F ) 

S ( I ) = S 1 

IF( I.NE.O) GO TO 457 

REWIND 36 

JJ=0 

C L 0 S E ( U N I T = 3 6 ) 

GO TO 12 
C L 0 S E ( U N I T = 3 1 ) 

CL0SE(UHIT=32) 
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C L 0 S E ( UNIT=35) 

STOP 

END 

£ * * * * * * * & * * * * * * * * * * * * * * * * * * * * * * * * * * & * * & & * & * * * * * * * £ £ * * * * * * * 

C LDF LO . F A : A LOAD FLOW PROGRAM 

********************************************************* 



DIMENSION 

E( 33) , Y1 ( 3 3 , 33) , PSCH( 33) , QSCH( 33) , PCAL( 33) , 



1QCAL( 33) , PDIFF( 33) ,QDIFF( 33) , IC( 33) , S( 33) , JAC(68 , 68) 

) 

2PQVEC( 33) , EFDV(6A ) , QQVEC( 33) , 

3YS( 33 , 33) , Y(33, 33) , CRFLO( 33 , 33) , U(68) , V(68) , 
AEDIFF(33) 

INTEGER P,Q,M,U,V,PP 

COMPLEX E,IC,Y,Z,YS,Yl,PFLO,PLOS,YSH,EDIFF, 
IS, CRFLO, PFLOl ,ZG1 , ZNI , ZNJ 
EQUIVALENCE (CRFLO, Yl) 

REAL JAC 

DIMENSION E I DV ( 6 A ) 



OPEH(UNIT=30,DEVICE= , DSK' , FILE =/ TEST. DAT' , DISPOSE='SA 
VE' ) 

OPEN( UNIT- 3 3 , DEVICE- ' DSK' , FILE= ' DSK. XFR' ) 



* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 5V * * 



CREAD VALUES OF IMPEDANCES YSIIUNT ADMITTANCES 



C 888888***** IF SHUNT ADMITTANCES ARENOT IGNORED 

REMOVE 

* * * T H E C IN THE LINE YS(Q,P) *********** *************** 

C***************** ************************** ********** 

READ( 30,2) IB, El 

2 FORMAT ( I , F) 

E( 1 )=CMPLX(E 1 ,0.0) 

DO 10 1 = 1 , 1000 

READ(30 , 5)NN, P, Q, J J , Z , ZG1 , ZNI , ZNJ 

3 FORMAT( A I , 8 F) 

IF(P.EQ.O) GO TO 15 

Y S ( P , Q ) = Y S H 

C YS(Q,P)=YSH 

IF(CABS(Y1(P,Q)).NE.0)Z=1.0/(Y1(P,Q)+1/Z) 

Y1(P,Q)=1 .0/Z 
Y1(Q,P)=Y1(P,Q) 

PRINT 6,P,Q,Y1(P,Q),YS(P,Q) 

6 FORMATC 21 , AF) 

10 CONTINUE 

********************************************************** 
CCOMPILE TERMS FOR YBUS 

Q********************* ******************************** 

15 DO 5 0 P = 1 , 1 B 

Y ( P , P ) = 0 
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DO 40 0=1, IB 
IF(P.EQ.Q) GO TO 40 
Y(P,P)=Y(P,P)+Y1(P,Q)+YS(P,Q) 
Y(P,Q)=-Y1(P,0) 

40 CONTINUE 

50 CONTINUE 

£*********************************************** 

C PRINT OUT VALUES OF YBUS 

Q***************** ******************************** 



PRINT 51 

51 FORMAT ( 'O' , 10X, ' BUS IMPEDANCE MATRIX , YBUS ' ) 

DO 58 1=1 , IB 

PRINT 54 , (Y( I , J) , J=1 , IB) 

54 FORMAT( '0' , 100F) 

58 CONTINUE 



£************** * * * * * * * * * * ************ ************** 



C READ VALUES OF POWER AT BUSSES 
£************************************************** 



DO 60 P=1 , IB 
READ( 30 , 59 ) I,PS,QS 
I F ( I . EQ . 0 ) GO TO 62 
IF(I.EQ.1)G0 TO 60 
PSCH( I )=-PS 
QS CH( I ) = -QS 

59 FORMAT ( I , 2F ) 

60 CONTINUE 



£************************************************** 



C ASSUME VALUES OF VOLTAGE AT BUSSES 

C************ ****************************** ******** 



62 DO 70 P=2 , I B 

E( P)=( 1 .0,0.0) 

70 CONTINUE 

K= 0 

DO 80 1=1 , IB 

TYPE 7 5 , I , E ( I ) 

75 F0RMAT(I,2F) 

80 CONTINUE 

122 PRINT 1 2 3 , K 

123 FORMAT (' 0 85X ITERATION = ',12) 

Q*********************************** ************* 



C CALCULATES REAL REACTIVE BUSS POWERS 
c* ************************************* ********** 



DO 150 P= 2 , I B 
P S U M = 0 
Q S U M = 0 

DO 140 Q=1 , IB 



P 1 = REAL( E( P ) )*( REAL(E( Q) )*REAL(Y(P , Q ) ) - A I M AG ( E ( Q ) ) * A I 
MAG ( Y 



1(P»Q))) + AIMAG(F.(P))*(AIMAG(E(Q))*REAL(Y(P,Q) ) + rcal( E 
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( Q ) ) * 

2 A I M AG (Y(P,Q))) 

Q1=AIMAG(E(P))*(REAL(E(Q))*REAL(Y(P,Q))-AIMAG(E(Q))*A 
I MAG ( Y 

1 (P ,Q) ) )-REAL( E( P) ) *( AIMAG(E(Q) ) * RE AL ( Y ( P , Q ) ) + RE AL ( E ( 
Q))*AI 

2MAG(Y(P,Q))) 

P S UM = P S UM+P 1 
QS UM = Q S UM+Q 1 
140 CONTINUE 

P C A L ( P ) = PSUM 
QCAL( P)=QSUM 

PRINT 145 , PCAL(P) ,QCAL(P) 

145 FORMAT( 2 F ) 

150 CONTINUE 

Q * * * * * *********************************************** 

C PRINT OUT VALUES OF CALCULATED BUSS POWERS 

Q* ****************************************** ********* 

PRINT 155, K 

155 FORMAT( ' 0 ' , ' REAL AND REACTIVE BUSS POWERS K=',I2) 

DO 158 1=2, IB 
PRINT 157 ,PCAL(I) ,QCAL(I) 

157 FORMAT ('0',2F> 

158 CONTINUE 
M= 0 

DO 160 P = 2 , IB 
EP=0 .01 

P D I F F ( P ) = P S C H ( P)-PCAL(P) 

QDIFF(P)=QSCH(P)-QCAL(P) 

TDIFF = ABS( PDIFF( P) ) + ABS( QDIFF( P) ) 

IF(TDIFF.LT.EP) CO TO 160 
M= 1 

160 CONTINUE 

C************** *************************** ****** 

C PRINT OUT VALUES OF DIFFERENTIAL BUSS POWERS 

C* ************************** ********************** 

PRINT 1 6 5 , K 

165 FORMAT( ' 0 DIFFERENTIAL PQ K=',3X,I2) 

DO 168 1=2,13 

PRINT 167 , PDIFF( I) ,QUIFF( I) 

167 FORMAT ('0',2F) 

168 CONTINUE 
IF(M.EQ.O) GO TO 400 
DO 170 I =2, IB 

S( I ) = CMPLX( PCAL( I ) , QCAL( I ) ) 

170 CONTINUE 

c* ********************************** ********* 

C CALCULATE BUSS CURRENTS I(P) EXCEPT AT SLACK BUS 
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C* *********************************** ******** 

1)0 175 1=2 , IB 

IC( I)=CONJG(S( I ) )/CONJG( E( I) ) 

PRINT 1 7 3 , I C( I ) 

173 FORMAT ( F ) 

175 CONTINUE 

c************** ********************** ****** 

C CALCULATE JACOBIAN USING PARTITIONS ' J 1 , J 2 , J 3 , J 4 ' 
c* *********** ******************** ********** 

C JACOBIAN 'Jl' 

C*************** 

IBDJ=2*IB-2 
I B J = I B- 1 
DO 190 M= 1 , I B J 
DO 180 N =1 , I B J 
P = M+1 
Q-N+l 

JAC(M,N)=(REAL(E(P))*REAL(Y(P , Q ) ) ) + ( A I M AG ( E ( P ) ) * AIM AG 
(Y(P,Q))) 

IF(P.NE.Q) GO TO 177 
JAC(M,N)=JAC(M, N)+REAL( IC( P) ) 

177 CONTINUE 

C PRINT 178,JAC(M,N) ,E(P) ,Y(P,Q) 

178 FORMAT ( 5 F ) 

180 CONTINUE 

190 CONTINUE 

*********************************************************** 
C JACOBIAN ' J 2 ' 

********************************************************** 
DO 210 M= 1 , IBJ 
DO 200 N=IB, IBDJ 
P = M+1 
Q=N-IB+2 



JAC(M,N)=-REAL(E(P) ) * A IM AG ( Y ( P , Q ) )+ AIM AG ( E ( P ) ) * RE AL ( Y 
( P > Q) ) 



IF(P.NE.Q) GO TO 195 
J A C ( M , N ) = J A C ( M , N ) + A I M A G ( I C ( P ) ) 
195 CONTINUE 

200 CONTINUE 

210 CONTINUE 



********************************************************* 



C JACOBIAN ' J 3 ' 

******************************************************* 



DO 230 M = I B , IBDJ 
DO 2 20 N = 1 , I B J 
P=M- I B+2 
Q = N+1 



JAC(M , N)=-REAL(E( P) )*AIMAG( Y( P , Q) )+ AIM AG ( E( P) )*R EAL( Y 
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( p ♦ Q ) ) 

IF(P.NE.Q) GO TO 215 

JAC( M , N ) = J A C ( M , N ) - A I M AG ( IC( P) ) 

215 CONTINUE 

220 CONTINUE 

230 CONTINUE 

ft ft ftftftftftftftftftft ft ftftftftftft ft ftftftftftftftftftft ftftftftftftftftftft ft ftftftftftftftftftftftftftftftft 

C JACOBIAN ' J 4 ' 

ftftftftftftftftftft ft ft ftftftftftftftftftftftftftftftftftftftftftftft ft ftftftftftftftftftftftftftftftftftft 

DO 250 M = I B , IBDJ 
DO 240 N = I B , IBDJ 
P=M- I B+2 
Q= N— I B + 2 

J A C ( M , N ) = — A I M AG ( E ( P ) ) * AIM AG ( Y ( P , Q ) )-REAL(E(P) )*REAL(Y 



( P > Q ) ) 


IF(P.NE.Q) GO TO 235 

JAC( M, N)= JAC( M, N) + REAL( IC( P) ) 


2 35 


CONTINUE 


240 


CONTINUE 


250 


CONTINUE 
DO 334 M= I , IBDJ 


C 


PRINT 333, (JAC(M,N) ,N=1 , IBDJ) 


333. 


FORMAT ( I 8 F ) 


334 


CONTINUE 



ftft ftftftft ft ftftftftftftftftftftftftftftftftftftftftftftftft ft ftft ft ftftftftftftftftftftftftftftftftftftftftftft ft ftft 
ft ft ftftftftftftftftftftftftftftftftftftftftftftftftftft ft ftftftftftftftftftftftftftftftftftftftftft ft ft ftftftft ftftftft 

C COMPUTE INVERSE OF JACOBIAN 
CALL MINV( JAC, 68 , IBDJ) 

ft ft ftftftftftftftftftftftftftftftftftft ft ftft ftftftftftftft ft ftftftftftftftftftftftftftftftftftft ftft ftftftftftft 

C PRINT 344 ,( JAC( M,N) ,N=1 , IBDJ) 

344 FORMAT( 1 3F) 

345 CONTINUE 

ftftftftftftftftftftftftftftftftftftftftftftftftftftft ft ft ftftftft ft ftftftftftftftftftftftftftftftftftftftftftftftft 

C COMPUTE DIFERENCE IN VOLTAGE USING JACOBIAN INVERSE 'JAC' 

ftft ft ftftftftftft ftftftft ft ftftftftftftftftftftftft ft ftftftftft ft ftft ft ft ftftftftftftftftftftftftftftftftftftft ft 

DO 346 1=2, IB 
M= I - 1 

PQVEC(M)=PDIFF( I ) 

Q Q V E C ( M)=QDIFF( I ) 

346 CONTINUE 
J = I B— 1 

DO 351 M = 1 , J 
K F D V D = 0 
E F D Q D = 0 
EIDD=0 
EIQV=0 

DO 34 9 N = 1 , J 

EFDVD=EFDVD+JAC( M , H ) * P Q V E C ( N ) 

349 CONTINUE 

DO 350 N= IB, IBDJ 



102 



NQ=N-J 

EFDQD=EFDQD+JAC( M, N)*QQVEC( NQ) 

350 CONTINUE 
EFDV(M)=EFDQD+EFDVD 

351 CONTINUE 

DO 360 M = I B , IBDJ 
E 1 DD=0 
E 1 QV= 0 

DO 3 58 N = 1 , J 

E1DD=EIDD+JAC( M , N) *PQVEC( N) 

358 CONTINUE 

DO 359 N=IB,IBDJ 
NQ=N-J 

EIQV=EIQV+JAC(M,N)*QQVEC( NQ) 

359 CONTINUE 
EIDV(M)=EIQV+EIDD 

360 CONTINUE 

* ******** * ************************ ********** ************* 

C COMPUTE NEW BUS VOLTAGES USING ' E F DV ' EFDIF FERENTI AL 
VECTOR 

******************************************************** 

DO 370 P = 2 , 1 B 
M=P- 1 
N = M+J 

IBZ=2*( IB-1 ) 

IF(N.GT.IBZ) GO TO 370 

EDIFF( P)=CMPLX(EFDV(M) , EIDV(N) ) 

PRINT 368 , P , EDIFF ( P ) 

368 FORMAT ( I , 2F) 

E ( P) = E(P)-rEDIFF(P) 

370 CONTINUE 

IF( K.EQ. 1 0 ) GO TO 390 
K=K+ 1 
GO TO 122 
390 TYPE 395 

395 FORMAT ( AN ERROR IS MADE IN DATA, LOAD FLOW DID NOT 

CONVERGE' ) 

A 00 CONTINUE 

500 CLOSE ( UNIT=30) 

********************************************************* 

C COMPUTE POWER FLOWS BUS P TO Q 

C COMPUTE CURRENT FLOWS AND POWER LOSSES 

C COMPUTE SLACK BUS POWER 

*********************** * *************************** ****** 

OPEN( UNIT=3 1 , FILE= ' DSK' ) 

S( 1 )=0 

503 FORMAT ( I ) 

DO 520 P= 1 , I B 
DO 515 Q=1 , IB 
IF(P.EQ.Q) GO TO 515 
IF(Y(P,Q).EQ.(0. 0,0.0)) GO TO 515 
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PP = P 

*Vc 'k'k-k*X-k*;’k‘k'k'kX'k , kJc'k'k-k'k*'kxi<&'k'k'kX>\*-k&'k&-k-k'k%'kX'k**&'k'k&'k i i( 

C IF YOU ARE NOT NEGLECT SHUNT ADMITT REMOVE THE NEXT 

C 

&£*&***£*** i; &&****&** >V *&**********£ ***** vc & ‘/j^^x^**** 

PFL0=C0NJG(E(P))*((E(P)-E(Q))*Y1(P,Q)) 

C l+( CONJG( E( P) ) *E( P)*YS( P , Q) ) 

PFL01=C0NJC(E(Q))*((E(Q)-E(P))*Y1(P,Q)) 

C 1 + ( COMJG(E( Q) )*E(Q)*YS( P , Q) ) 

CRFLO( P , Q)=PFLO/ CONJC(E(P) ) 

PLOS = PFLO+P F LO 1 

********************************************************* 

C TAKE THE CONJUGATE OF THE POWER FLOWS FOR FORMAT 

P-JQ 

********************************************************* 
PLO S = CO N JG ( PLOS) 

PFLO=CONJG( PFLO ) 

IF(P.NE.l) GO TO 505 
S( 1 ) = P F LO+ S ( 1 ) 

505 WRITE( 31 , 509 )P , Q, PFLO, CRFLO( P , Q) , PLOS 

WRITE(31,510)PP,E(PP) 

509 F0RMAT(6X, ' POWER FLOW 

BUS' , 12 , ' TO' , 12 , ' 1 S' , 2F, 2X, ' CURRENT 

1FL0W IS' , 2F, ' POWER LOSS IS',2F) 

510 F0RMAT(5X, ' BUS' ,12, 'VOLTAGE 1S',2F) 

515 CONTINUE 

520 CONTINUE 

WRITE ( 3 1 , 5 AO) S ( 1 ) 

540 F0RMAT(6X, ' POWER AT SLACK BUS 1 IS',2F) 

DO 550 M= 1 , IB 

DO 545 N = 1 ,13 

I F ( CABS(CRFLO(M,N) ) .GT. 0.000001) 

W R 1 T E ( 3 3 , 5 4 3 ) M , N , 

1CRFL0( M, N) 

543 FORMAT ( 2 1 , 2 F ) 

545 CONTINUE 

550 CONTINUE 

WR1TE( 3 3 , 555 ) 

555 FORMAT( ' 0 0 0.0 0.0') 

CLOSE( UN1T=33) 

DO 600 M=1 , IB 
DO 5 80 N = 1 , I B 

S( 2 )=E( M)*CONJG(CRFLO(M, N) ) 

PFAD=CABS( S( 2 ) ) 

IF(PFAD.EQ.O) GO TO 580 
PFAC— ( REAL( S( 2 ) ) )/ ( PFAD) 

WR I TE ( 3 1 , 5 7 0 ) M , N, PFAC 

570 FORMAT( 4X, ' POWER FACTOR BU S ' , I 3 , ' TO ' , I 3 , ' 1 S ' , F ) 

580 CONTINUE 

600 CONTINUE 

CLOSE ( U N I T = 3 1 ) 



c 



c 

c 

c 



10 

15 



20 



25 



3 0 
35 

38 



AO 

45 

50 

55 



6 0 
62 
65 



70 

75 

8 0 

100 



104 



STOP 

END 

TITLE: MI NV;SUBR0UT1NE TO INVERT A REAL SQUARE MATRIX 
SUBROUTINE MINV( A, MDIM , N) 

DIMENSION A( MDIM, K) 

INTEGER L( 100) ,M( 100) 

CONVERTED FROM SSP ROUTNE MINV 

DO 80 K= 1 , H 
L(K) = K 
M ( K ) = K 
BIGA=A( K, K) 

DO 20 J = K , N 
DO 2 0 I = K , K 

IF( ABS( BIGA)-ABS( A( I , J ) ) ) 15,20,20 
BIGA=A( I , J ) 

L(K)=I 
M( K)=J 
CONTINUE 
J = L ( K ) 

IF(J-K) 35,35,25 
DO 30 1 = 1 , N 
HOLD=-A( K, I) 

A(K,I)=A( J, I) 

A( J , I ) = H 0 L D 
I = M ( K) 

IF(I-K) 45,45,38 
DO 4 0 J = 1 , N 
HOLD=- A( J , K) 

A(J , K ) = A( J , I) 

A( J, I ) = H 0 L D 

DO 55 1=1 , N 

IF(I-K) 50,55,50 

A( I , K)=A( I , K) / (-BIGA) 

CONTINUE 
DO 65 1=1 , N 
HOLD=A( I , K) 

DO 65 J = 1 , N 
IF(I-K) 60,65,60 
IF(J-K) 62,65,62 
A ( I, J)=HOLD*A(K,J)+A( I ,J) 

CONTINUE 
DO 75 J=1 ,N 
IF(J-K) 70,75,70 
A(K, J)=A( K, J) / BIGA 
CONTINUE 

A(K, K) = ( 1 . 00 00000) / BIGA 
CONTINUE 
K = N 
K=K- 1 
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105 


I F ( K ) 150,150,105 

i=l(k) 


108 


I F ( I-K) 120,120, 108 
DO 110 J = 1 , N 


1 10 


HOLD=A( J , K) 

A( J,K)=-A( J, I) 
A( J, I ) = HOLD 


120 


J = M ( K ) 


125 


IF(J-K) 100,100,125 
DO 130 1=1, N 


1 30 


HOLD=A( K , I ) 
A(K, I ) =-A( J , I) 
A( J , I ) = H 0 L D 


150 


GO TO 100 
RETURN 




END 



****************************************************** 

C TITLE : ZBUS • F 4 ;ACCEPTS DATA FROM TEST. DAT AND DSK.XFR 

FROM 
C 
C 

FAULT 
C 

CURRENTS . 

********************************************************* 
***** ******* ***************** **************************** 

C READ DATA FROM TEST. DAT 

********************************************************* 

DIMENSION 

ZD( 300 , 3) , ZM( 300 , 2 ) , RID( 30 , 30, 3) , ZG( 30) , ZN( 30) , 

110(200) , RI ( 3 0 ) ,ZB(3 0, 30) ,Z BO (30,30) , ZMU( 5 , 5) , IDD( 5) , 

ZL( 30) , 

2ZBR(30),ZG0(30),ZL0(30),EB(30) 

COMPLEX ZD,ZM,ZG,T1,ZN,ZN1,ZNI,ZNJ,ZG1,Z0, 

1ZB, ZBO, ZMU, ZBB, ZL, ZLL, ZBR, ZGO , ZLO , ZLLO , EB 
NN B= 1 

EQUIVALENCE (ZBR, ZLO , EB) 

OPEN (UNIT =31, FILE=' TEST. DAT ' ) 

C OPEN( UNIT=3 2, FILE=' DSK.XFR ' ) 

RE AD ( 3 1 , 10)IB,EB1 
10 F ORMAT ( I , F ) 

20 READ( 31 , 30) NB , I, J , JJ , T1 , ZG1 , ZNI, ZNJ 

30 FORMAT( 4 1 , 8F) 

IF(NB.EQ.O) GO TO AO 
DN = NB 
N N = N B 

Z D ( N N , 1 ) = CMPLX(1)N, 0 .0) 

Z 1) ( N N , 2 ) = T 1 
ZD( NN , 3 ) = Z G 1 + 3 * Z N ( I ) 

Z K ( I ) = Z N I 



LOAD FLOW PROGRAM LDFLO.FA AND FORMS ZBUS BY 
ALGORITHM THIS WILL BE USED TO COMPUTE ALL 

CURRENTS . LOAD FLOW COMPUTED MAXIMUM LOAD 
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ZN( J ) = ZNJ 
R I D ( I , J , J J ) = N N 
IF(NB.NE.O)GO TO 20 
40 DO 45 1=1 , IB 

C READ (32,43)T1 

43 F0RMAT(2F) 

43 CONTINUE 

DO 48 1=1 , IB 
C READ(32,46)J,T1 

46 FORM AT ( I , 2F ) 

EB ( J ) = T 1 

48 CONTINUE 

49 READ ( 31 , 50 ) I , T 1 , XD 1 

50 FORM AT( 1 , 3 F ) 

IF( I . EQ. 0) GO TO 60 
IF( CABS(T1 ) . EQ. 0) GO TO 49 
C T1 = CON J G ( T 1 ) 

C ZG(I)=3*(CABS(EB(I))**2)/T1 

1 F ( X D 1 . NE.0)ZG(I)=CMPLX(0.0,XD1 ) 
ZGO ( I ) =ZG( I ) 

IF(I.NE.O) GO TO 49 



********** * * * * * * ******* * *************** * * * * * * * * * * * * * * 



c 

* * * * * * * * 



ENTER THE ZERO SEQ IMF BUS TO REF( GROUND) 
************************************************ 



60 TYPE 62 

62 FORM AT ( ENTER THE SYNC ZERO SEQUENCE IMPUDENCE FOR 

THE LOADS 

1 AT EACH '/' BUS, 3*ZN, MUST BE ADDED, IMPED IN PU..,') 

DO 65 1=1 , IB 

ACCEPT 6 4 , 1 1 , T 1 

IF( II . EQ.O)GO TO 66 

ZGO ( 1 1 ) = T1 

64 FORMAT( I , 2F) 

65 CONTINUE 

66 DO 68 1=1 , IB 
IF(CABS(ZGO(I)).EQ.O)ZGO(I)=ZG(I) 

68 CONTINUE 

****** * * * * ********************************************** 



C ENTER THE MUTUALLY COUPLED ELEMENTS 

******************************************************** 

CLOSE( U K I T = 3 1 ) 

TYPE 70 

70 FORMAT ( ' ENTER THE MUTUALLY COUPLED ELEMENTS IN PER 

UNIT' / ' 

1 IN T1IE FOLLOWING FORMAT, TEST. DAT HAS THE ELEMENT 

S. . , '/' 

2 (ELEMENT )SP T 0 ( E L E M E N T )SP(0 SEQ VALUE) ') 

75 ACCEPT 80,EL,ELL,Z0 

80 FORMAT(4F) 

IF(EL. EQ.0.0) GO TO 100 
11 = 11+1 
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Z M ( II, 1 ) = CMPLX( EL, ELL) 

ZM( II, 2 ) = Z 0 
GO TO 75 

*&***&&#*£**£**£***£****# s'c **&&** & v'c***#*********}^*****:!? 

C ALL ELEMENT DATA PREPARED PERFORM POSITIVE AND 

NEGATIVE 

C SEQUENCE COMPUTATIONS TO FORM ZBUS 

***********&************************************&******* 
******************************* * ********************* 



100 



120 

no 

115 



no 

1 1 5 



1 15 



120 



DO 500 I EL= 1 , NN 
BRA = 0 
ID( I EL ) = 1 
DO 120 1=1 , 30 

DO 120 J=1 , 30 

DO 1 20 JJ=1 , 3 

IF(RID(I, J, JJ) .NE.IEL) GO TO 



CONTINUE 



IF(RI(I) . EQ . 0 ) GO TO 

IF( RI( J) .EQ.O)GO TO 

BRA= 1 . 0 

IPP = MIN0( I , J ) 

I Q Q = M A X 0 ( I , J) 

GO TO 130 

IF(RI(J) . E Q . 0 ) G 0 TO 

RI( I)=l .0 
IPP = J 
IQQ=I 
GO TO 130 
R I ( J ) = 1 .0 
R I ( I ) = 1 . 0 
I P P = MI NO ( I , J) 
IQQ=MAX0 ( I , J ) 

GO TO 130 



***************************************************** 



C BRANCH OR LINK HAS BEEN DETERMINED BRA=0 : BRANCH 

********** * ************************** A ******************* 



130 DO 195 I L= 1 , I B 

ZBB=CMPLX(0. 0,0.0) 

181 IF( BRA.NE.O) GO TO 186 

I F ( IQQ.EQ.IL) GO TO 190 
Z BR ( IL) = ZB( IPP, IL) 

GO TO 190 

186 ZL( IL)=ZB( IPP, IL)-ZB( IQQ, IL) 

190 CONTINUE 

195 CONTINUE 

Z BB= CMP LX ( 0 .0,0.0) 

I F( BRA.NE.O) GO TO 250 
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DO 200 1 L = 1 , IB 

Z B ( IQQ, IL) =Z BH( I L) 

ZB( I L , 1QQ)=ZBR( IL) 

200 CONTINUE 

240 ZB(I0Q,IQQ)=ZB(IPP,IQ0)+ZD(IEL,2) 

Z B B = C M P L X ( 0 . 0 , 0 . 0 ) 

GO TO 300 

250 ZLL=ZL( 1PP)-ZL( IQQ)+ZD( 1EL, 2) 

ZBB=CMPLX(0 .0, 0.0) 

DO 290 1 P= 1 , IB 
DO 290 I Q= 1 , IB 

I F ( IP.GT. IQ)GO TO 290 



ZB(IP,IQ)=ZB( IP, IQ)— (ZL( I P ) * Z 
L(IQ) 

1 /ILL) 

IF(IP.EQ.IQ) GO TO 290 
Z B ( IQ, IP) = ZB( IP, IQ) 

290 CONTINUE 

300 DO 350 1=1,5 

DO 350 J = 1 , 5 

ZMU ( I , J ) = CHPLX ( 0 . 0 , 0 . 0 ) 

I DD ( I)=0 

350 CONTINUE 

DO 400 1=1 , IB 
DO 400 J = 1 , IB 

PR- 1 NT 380, I, J,ZB( I, J) 

380 FORMAT( 2 1 , 2 F ) 

400 CONTINUE 

NNB= 1 

500 CONTINUE 



-k'k'k-k'k'k'k'k'k'k It •k'k'k'k'k'k'k'k'k'k-k'k'k-k-k'k'it'k'k'k'k'k'k'k'k-k'k-k'k'k'k-k-k-k'k-k-k'k'k'k'k 



c 



ZBUS POS/NEC IS FORMED 



'k'k'k'k'k’k'k'k'k'k'k'k-k'k'k'k'k'b'k it 'k’kit'k'k'k'k'k'k'k'k'k'k'k'kit'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k 



DO 505 1=1 , NN 

I D ( I ) = 0 

505 CONTINUE 

DO 510 1=1, IB 

RI( 1 ) = 0 . 0 

510 CONTINUE 

DO 1000 I E L = 1 , N N 
BRA = 0 
I D ( 1 E L ) = 1 
N N B = 1 

DO 620 1=1,30 

DO 620 J=1 , 30 

DO 620 JJ= 1 , 3 

1F(RI1)( 1 , -J , J J) . NE. I EL) GO TO 

620 



6 1 0 



IF(RI(I). EQ. 0) GO TO 
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IF(RI(J) . E Q . 0 ) GO TO 

615 

BRA = 1 . 0 
I P P = M I N 0 ( I , J ) 

I QQ=MAX0 ( I , J) 

GO TO 625 

610 IF(RI(J).EQ.O) GO TO 

6 1 5 

RI ( I ) = 1 .0 
IPP = J 
IQQ=I 
GO TO 625 

615 R I ( I ) = 1 .0 

R I ( J ) = 1 . 0 
IPP=I 
IQQ = J 
GO TO 625 

620 CONTINUE 

&&**&&&&&&&*&*&& & ******* & >*c ****** & * ****£*****&&****** >Wc * 

C SEARCH FOR MUTUAL COUPLING IN ZERO SEQUENCE 

625 IDD( 1 )=IEL 

DO 650 1=1 , NN 
DM=I EL 

IF(REAL(ZM(I,1)) . EQ . DM ) MUT= 1 
I F ( A I M AG ( Z M ( I , 1)) . EQ. DM)MUT=1 
IF(MUT. ME. 1 )GO TO 650 
M U T = 0 

II 1=IFIX( REAL(ZM( I , 1))) 

1 1 2= I F IX ( AIMAG( ZM( I , 1))) 

I F ( I D ( 1 1 1 ) .EQ.O)GO TO 650 

IF( ID( I 12) . EQ. 0)GO TO 650 

IF( 1 1 1 .EQ.IEL) GO TO 630 

113=11 1 

111=112 

112=113 

630 I D D ( 1 ) = 1 1 1 

NNB=NNB+1 
I DD( NNB) = 1 12 
DM=I 12 

DO 645 J 1 = 1 , N H 

I F( REAL( ZM( J 1 , 1 ) ) . EQ. DM) MUT=1 

I F ( A IM AG ( Z M( J 1 , 1 ) ) . EQ . DM ) MUT = 

1 

I F ( MUT . NE . 1 ) GO TO 64 5 
M U T = 0 

I I I 1 = I F 1 X ( R E A L ( Z M ( J 1,1))) 

I I I 2 = 1 F 1 X ( A I M AG ( Z M ( J 1 , 1 ) ) ) 

I F ( I D ( I I 1 1 ).EQ.O)GO TO 64 5 
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IF(1D(1I12).EQ.0)G0 TO 645 
IF( I III .EQ. DM) GO TO 64 3 
1113=1111 
1111=1112 
1112=1111 

643 IF( III2.EQ. Ill ) GO TO 645 

NNB=NNB+1 

IDD(NNB)=III2 

645 CONTINUE 

650 CONTINUE 

DO 660 1=1 , NNB 

Z MU ( I , I ) = ZD( IDD( I ) , 3) 

DO 660 J=1 , NNB 

DD1=IDD( I ) 

DD2=IDD(J ) 

DO 660 KK= 1 , NN 



IF(ZM( KK, 1 ) . EQ. CMPLX( 
DD1 , DD 2 ) ) 

1 M U T = 1 



IF(ZM( KK, 1 ) . EQ. CMPLX( DD2 , DD 1 ) 
) M U T = 1 



660 



CONTINUE 



IF( MUT. NE. 1 )G0 TO 660 
MUT = 0 

ZMU ( I , J) = ZM(KK, 2) 
ZMU( J , 1 )=ZM( KK, 2 ) 



********************************************************* 



C ZMU IS FORMED 

* ***** * ********** *********** * ************************** 



DO 665 1=1 , NNB 
DO 665 J=1 , NNB 
PRINT 663, I, J, ZMU (I, J) 

663 FORMAT (21, 2 F ) 

665 CONTINUE 

********************************************************* 



CALL CMXINV( ZMU , 5 , NNB) 



******************************************************** 



DO 663 1=1, NNB 

PRINT 667 , ( ZMU( I , J ) , J=1 , NNB) 

667 FORMAT ( 1 OF ) 

663 CONTINUE 

DO 695 1 L = 1 , IB 
IF(NNB.EQ.l) CO TO 631 
ZBB=CMPLX( 0 .0,0.0) 

DO 680 J 4 = 1 , N N B 
J = J 4+1 
I = I DD ( J ) 

IF(I.EQ.O) GO TO 630 
DO 670 11=1,30 



Ill 



13=1,3 
1 F ( R I 

0(11,12,13 

1). HE. I) GO TO 670 



DO 670 12=1,30 
DO 67 0 



IP = I 1 



I Q= 1 2 

ZBB=(ZMU(1 ,J)/ZMU( 1 , 1) )*(ZB0(IP, IL)-ZB0(IQ, IL) )+ZBB 
670 CONTINUE 

680 CONTINUE 

681 IF(BRA.NE.O) GO TO 686 

IF (IQQ.EQ.IL) GO TO 690 
ZBR(IL)=ZBO(IPP,IL)+ZBB 
GO TO 690 

686 ZL( IL)=ZBB+ZB0( IPP, IL)-ZBO(IQQ,IL) 

690 CONTINUE 

695 CONTINUE 

Z B B= CMPLX ( 0 .0,0.0) 

IF(BKA.NE.O) GO TO 750 
DO 700 I L= 1 , IB 

ZBO ( I QQ , IL)=ZBR(IL) 

Z BO ( 1 L , IQQ) = ZBR( IL) 

700 CONTINUE 

IF( NNB . EQ . 1 ) GO TO 740 
DO 730 J=2 , NNB 

I=IDD(J ) 

DO 7 20 11 = 1 ,30 

DO 720 12=1,30 

DO 720 13=1 , 3 

IF(RID(I1 ,12, 13) .NE.I) GO TO 720 

IP=I 1 
IQ= 1 2 

720 CONTINUE 



ZBB=(ZMU( 1 , J )*(ZB0( IP, IQO)-ZBO( IQ, IQQ 

))/ 

1 Z M U ( 1 , 1 ) )+ZBB 
730 CONTINUE 

740 ZBO (IQQ, IQQ)=ZB0( I PP , IQQ)+ZBB+1 . 0/ZMU( 1 , 1 ) 

ZBB = CMPLX(0. 0,0.0) 

GO TO 800 

750 IF(NNB.EQ.l) GO TO 781 

DO 7 80 J = 2 , NNB 

I=IDD(J) 

DO 770 11=1,30 

DO 770 12=1 , 30 

DO 7 7 0 13= 1 , 3 
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1 F ( R 1 D ( 1 1 , 

112 , 13) . NE. 1 )C0 TO 7 70 

IP = I 1 
I Q = 1 2 

770 CONTINUE 

Z BB = ( Zi-lU ( I , J)*(ZL( I P)-ZL( IQ) )/ZMU( 1 , l 
) )+ZBB 

780 CONTINUE 

781 ZLL=ZL(IPP)-ZL(IQQ)+ZBB+(1.0/ZMU(1,1)) 

ZBB=CMPLX(0. 0,0.0) 

DO 790 1P=1 , IB 
DO 790 I Q= 1 , IB 

IF( IP.GT. 1Q)G0 TO 790 

ZBO(IP,IQ)=ZBO(IP,IQ)— (ZL(IP) 

*ZL( IQ) 

1 /ZLL) 

IF(IP.EQ.IQ) GO TO 790 
ZBO ( I Q , I P ) = Z BO ( 1 P , IQ) 

790 CONTINUE 

800 DO 850 1=1,5 

DO 850 J=1 , 5 

ZI1U(1,J)=CMPLX(0. 0,0.0) 

1DD( 1 )=0 

850 CONTINUE 

DO 900 1=1 , IB 
DO 900 J=1 , IB 

PRINT 880, I, J, ZBO(I , J) 

080 FORMAT ( 2 I , 2 F ) 

900 CONTINUE 

N N B = 1 

1000 CONTINUE 

************************************************************* 

****** 

C Z BUS ZERO IS FORMED 

**************************************************** 

OPEN(UNIT=33 , FILE='ZBUS . FOR' ) 

DO 1010 1=1, IB 
DO 1010 J=1 , IB 

W R I T E ( 3 3 , 1005)1, J, ZBO (I, J),ZB(I,J),ZB(I,J) 

1005 FORMAT( 21 , 6F) 

PRINT 1007 

1 007 FORMAT ( 3 X , 'BUS TO BUS ',3X,'ZER0 SEQ ',5X,'P0S 

SEQ' ) 

PRINT 1008,1, J,ZBO(I, J) ,ZB(I, J) 

1008 FORMAT ( 2 1 6 , 4 F ) 

1010 CONTINUE 

CLOSE(UNIT=32) 

**************************************************** 

C NOW TO ADD ALL LOADS. ALL LOADS WILL BE LINKS WITH 
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NO 

C MUTUAL COUPLING. 

*********************************************** * ****** 
DO 1200 1 = 1 , I B 

IF(CABS(ZG(I)).EQ.O)GO TO 1200 
DO 1100 I L= 1 , I B 

Z L ( IL)=-ZB( I , IL) 

ZL0(IL)=-ZB0(I , IL) 

1100 CONTINUE 

ZLL = -ZL( I ) + Z G ( I ) 

ZLL0=-ZL0( I )+ZG0( I ) 

DO 1150 J = 1 , IB 

DO 1150 K=1 , IB 

I F ( J . GT . K ) GO TO 1150 

ZB(J,K)=ZB(J,K)-ZL(J)*ZL(K)/Z 

LL 

ZB(K,J)=ZB(J , K) 

ZB0(J,K)=ZB0(J,K)-ZL0(J)*ZL0( 

K ) / Z LLO 

ZBO (K , J )=ZB0( J , K) 

1150 CONTINUE 

1200 CONTINUE 

OPEN(UNIT=34,FILE=' ZBUS.DAT') 

DO 1 300 1=1 , IB 

DO 1300 J = 1 , I B 

WRITE(34,1250)I,J,ZB(I,J),ZB0(I,J) 

1 2 50 FORMAT ( 2 1 , 4 F ) 

1300 CONTINUE 

CLOSE(UNIT=34) 

CLOSE(UNIT=33) 

OPEN(UNIT=36, FILE='MUT. CO' ) 

DO 1400 1=1, NN 

EL 1 = REAL ( ZM( 1,1)) 

IF(EL1.EQ.0.0)G0 TO 1322 
EL2=AIMAG(ZM( 1,1)) 

IEL1=IFIX(EL1 ) 

IEL2=IFIX(EL2) 

DO 1320 11=1 , IB 
DO 1320 J 1 = 1 , I B 
DO 1320 K 1 = 1 , I B 

IF(RID( I 1 , J 1 , K1 ) . EQ. ELI )GO TO 13 10 
IF( RID( I 1 , J 1 , K1 ) . EQ. EL2)GO TO 1315 
GO TO 1320 

1310 I B 1 1 = 1 1 

I B I J = J 1 
GO TO 1320 

1315 I B J I = I 1 

I B J J = J 1 

1320 CONTINUE 
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1322 IF (ELI . EQ.O)GO TO 1400 

WRITE (36, 1350)EL1 , E L 2 , I B 1 1 , IBIJ, I B J I , I B J J , Z t-1 ( 

1,2) 

1350 FORMAT ( 2F,4I3, 2F) 

1400 CONTINUE 

W R I T E ( 3 6 , 14 50) 

1450 FORMAT( ' 0.0 0.0 0 0 0 0 0.0 0.0') 

CLOSE(UMIT=36) 

STOP 

END 

C TITLE : CMXINV; SUBROUTINE FOR INVERT COMPLEX MATRIX 

SUBROUTINE CMXINV ( A , MDIM , N ) 

COMPLEX A ( M D I M , N ) , BIGA, HOLD 
INTEGER L ( 100) ,M( 100) 

C 

C CONVERTED FROM SSP ROUT NE MINV 

C 

DO 80 K= 1 , N 
L ( K ) = K 
M( K) = K 
BIGA=A( K , K) 

DO 2 0 J = K , N 
DO 2 0 I = K , N 

10 I F ( CABS( BIGA)-CABS(A( I , J ) ) ) 15,20,20 

15 BIGA=A( I , J ) 

L ( K ) = I 
M( K) = J 

20 CONTINUE 

J=L(K) 

IF(J-IC) 35,35,25 
25 DO 30 1=1, N 

HOLD = - A( K , I ) 

A(K, I)=A( J, I) 

30 A(J,I)=HOLD 

35 I-M(K) 

IF(I-K) 45,45,38 
38 DO 40 J= 1 , N 

UOLD=-A( J , K) 

A(J ,K)=A(J , I) 

40 A( J , I ) = HOLD 

4 5 DO 5 5 1=1 ,N 

IF(I-K) 50,55,50 
50 A(I,K)=A( I , K ) / (-BIGA) 

55 CONTINUE 

DO 65 1=1 , N 
liOLD=A( I , K) 

DO 65 J= 1 , N 
IF(I-K) 60,65,60 
60 IF(J-K) 62,65,62 

62 A( I, J)=HOLD*A( K, J)+A( I , J) 



6 5 


CONTINUE 
00 7 5 J = 1 , n 


70 

75 


IF(J-k) 70,75,70 
A C K, J)=A( K, J)/BIGA 
continue 


80 


A(K,K) = (l. 00 00000 0 
CONTINUE ’ 




K=N 


100 


K = K- l 


105 


IF (K) 150,150,105 
1=L(K) 


108 


nn ( ?7n C) I 20 * 120 *1CS 
DO 110 J = 1 , N 

HOLD=A( J , K) 


1 10 


A ( J , K)=-A( J , I ) 
AC J , I ) = HOLD 


120 


J = MC K) 


125 


IF C J-K) 100, 100, 125 
D0 130 1=1 n 



1 30 
150 



HOLD=A(K, I) 

A( K, I ) =-A( J , I ) 
A ( J , I ) =HOLD 
GO TO 100 
RETURN 
END 



******************* a 

c c nm ir .LT.F4. Tm ; 

;;; 

1 Z R2 ^2N2 F DELE^FCUR^°fcu^^ * *ZNI, ZN J , 

-°S» •;* S -->’ ALCUR 

MAD( 34 ,5) ib"i ‘ : “ US -“«') 

5 E0RMAT(I,F) 

I E B= I B * I B 
D0 7 1 = 1 , IB 

EB(i)=cmp LX (E] ,o. 0 ) 

CONTINUE 

R E A D ( 3 4 , 10) , NI , I j y 

format ( 4 i ) ’ ,k 

IF ( NI . Eq . o ) Go TO 15 
IF(K.EQ. 2 )I=-i 

IF(k.EQ. 3)j=_j 
R I = I 
RJ = J 

I0ELE(NI )=CMPLX(RI , RJ) 



7 

9 

10 
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IF(NI.NE.O)CO TO 9 
15 DO 30 11 = 1 , IBB 

READ(35 , 25) I , J , ZB1 , Z B 0 1 
25 FORMAT ( 2 1 , 4 F ) 

Z B ( 1 , J ) = Z B 1 
Z BO ( I , J) = ZB01 
30 CONTINUE 

CL0SE(UN1T=35) 

OPEN( UNIT = 37 , FILE='MUT. CO' ) 

0PEH(UN1T=38, F I LE= ' F ALT . DAT ' ) 

0PEH(UH1T=39, FILE=' EFLT.DAT' ) 

OPE N( U N I T = 4 0 , FILE=' RELA. DAT' ) 

OPEN(UNIT=41 , FILE=' FCUR.DAT ' ) 

OPEN(UNIT=42, FILE= ' DSK. XFR' ) 

OPEN(UNIT=43, FI LE= ' LCUR . DAT' ) 

READ (40, 32) RMVAB , RKVBA, IRR 
3 2 FORMAT( 2 F , I ) 

RKVBA=RKVBA* 1 E3 

FIPU=RMVAB* 1 E 6 / ( SQRT( 3 .O) 5 - RKVBA) 

33 READ(40 , 34) IRR 

34 FORMAT(I) 

IF(IRR.NE.0)G0 TO 33 

36 READ ( 4 0 , 3 7 ) J , E 

37 FORM AT ( I , F ) 

IF(J.EQ.O)GO TO 3.3 
EIB(J )=E 

GO TO 36 

38 TYPE 35 

35 FORMAT ( ' ENTER THE FAULT 

IMPEDE MCE- REAL ( SP) IMAGINARY. . ) 

ACCEPT 40, ZF 
40 FORM AT( 2 F ) 

DO 500 I F L T = 1 , IB 

IF( CABS( ZB( IFLT, IFLT) ) . EQ. 0)GO TO 500 
DO 50 1=1 , IB 

IF( I . EQ. IKLT)GO TO 45 

E FALT ( 1 , 2, 1 ) = SQRT(3 . 0 ) * ( EB( I )- ( ( Z B( I , 

IFLT) *EB 

1 ( IFLT) ) / ( ZF+ZB( IFLT, IFLT) ) ) ) 

EFALTC I , 1 , 2)=-( SQRT( 3 . 0)*EB( IFLT)*ZB0 
(I, IFLT)/ 

1 (ZB0( IFLT, IF LT)+3*ZF+( 2*ZB( IFLT, IFLT) )) ) 

E FALT (I , 2 , 2 )=SQRT( 3 . 0)*EB( I )-( ( SQRT( 3 
. 0 ) * E B 

1 (IFLT)*ZB( I , IFLT) ) /( ZB0( IFLT, IFLT)+( 2*(ZB( IFLT, IFLT) 
)+3*ZF))) 
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E F ALT( 1 , 3 , 2) = EFALT( I , 2 , 2)-( SQRT( 3 .0)* 

EB( I ) ) 

GO TO 50 
45 

EFALT( IFLT, 2 , 1 )- SQRT( 3 . 0 ) * ( Z F* E B( IF LT 
) / ( Z F+ 

1 Z B ( IFLT, IFLT) ) ) 

E F ALT ( IFLT, 1 , 2)=-SQRT( 3 . 0 ) * EB ( I F L T) * Z 
B0( IFLT, 

1IFLT) / ( ZBO( IFLT, IFLT)+( 2 * Z B ( IFLT, IFLT) ) + 3 * Z F ) 

EFALT( IFLT, 2,2)=SQRT(3.0)*EB(IFLT)*(Z 
BO ( I F LT , 

1IFLT)+ZB( IFLT, IFLT)+3*ZF)/(ZB0(IFLT, IFLT)+(2*ZB( IFLT 
, IFLT) )+3* 

2ZF) 

EFALT( IFLT, 3 , 2) = -SQRT( 3 .0)*EB( IFLT)*Z 
B( IFLT, IFL 

IT) / (ZBO( IFLT, IFLT)+( 2 * Z B ( IFLT, IFLT) )+3*ZF) 

50 CONTINUE 

REWIND 3 4 
REWIND 37 

WRITE ( 39, 53) ( IFLT, I , E F ALT( 1,2,1) , E F A L T ( 1 , 1,2) 

IE F ALT (I, 2, 2) , E F ALT ( I, 3, 2), 1=1, IB) 

5 3 FORMAT ( 2 1 3 , 8 F ) 

READ (34, 55)11, E 
55 FORMAT ( I , F ) 

DO 90 K=1 , IBB 

READ( 34 , 6 0 ) N N , I , J , KK , ZR1 ,ZN1 ,ZNI,ZNJ 
60 FORMAT ( 4 I , 8 F ) 

IF(NN.EQ.O)CO TO 91 

FCUR( NN , 2 , 1 ) = ( EF ALT( 1,2, 1)-EFALT(J,2, 1 ) ) / Z R 1 

FCUR(NM, 1 , 2) = (EFALT( I, 1 , 2)-EFALT(J , 1 , 2))/ (ZN1 
+ 3 * 

1 ( Z N I + Z N J ) ) 

FCUR ( NN , 2, 2)=( EFALT( I , 2, 2)-EFALT(J , 2, 2) )/ZRl 
FCUR(NN,3, 2 ) = ( E F ALT( 1,3, 2 )-EFALT( J , 3 , 2))/ZRl 
DO 85 12=1 , IBB 

READ( 37 , 65)RN1 , RN2 , IBI I , IB I J , 

IBJI , IBJJ 

1 , ZN 2 

6 5 F 0 RM A T ( 2 F , 41 , 2 F ) 

N N 1 = I F I X ( R N 1 ) 

N N 2 = 1 F I X ( R N 2 ) 
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1F( NN 1 . EQ. MU) GO TO 70 
1 F ( N N 2 . E Q . N N ) GO TO 7 4 
I F ( N N 1 . E Q . 0 ) C 0 TO 8 6 
GO TO 85 



70 


FCUR(NN, 1 , 2) = ( ( EFALT( IBJI , 1 , 2 )- E F ALT ( IBJJ , 
1 (ZN2) )+FCUR(NN, 1 , 2) 

GO TO 85 


1,2))/ 


74 


FCUR(NN, 1, 2) = ( ( E F AL T ( I B 1 1 , 1 , 2 ) - E F ALT ( I B I J , 
1 ( Z N 2 ) ) + FCUR( NN, 1,2) 


1,2))/ 


85 


CONTINUE 




86 


REWIND 37 




90 


CONTINUE 




91 


REWIND 34 
DO 300 NN= 1 , IBB 





IF(CABS(FCUR(NN,2, 1)).EQ.O)CO TO 300 

WRITE (38, 250)1 FLT,NN, FCUR( NN , 2, 1) , FCU 
R(NN, 1,2) 

1 , FCUR(NN, 2, 2) , FCUR( NN, 3 , 2) 

250 FORMAT( 2 1 3 , 8 F ) 

300 CONTINUE 

DO 400 NN= 1 , IBB 

IF( CABS(FCUR(NN, 2, 1 ) ) . EQ. 0)GO TO 400 

FCUR1 = (RKVBA/EIB( I F IX ( AB S ( REAL ( I DELE ( 

NN) 

1)))))*FCUR(NN, 2, 1) 

FCUR2=(RKVBA/EIB( I F I X ( A B S ( RE AL ( I DE L E ( 

NN) 

1 ) ) ) ) )*(FCUR( NN, 1 , 2)+FCUR( NN, 2 , 2)+FCUR(NN, 3 , 2) ) 

FCUR1=FCUR1*FIPU 
FCUR2=FCUR2*FIPU 
IR 1= ( I F I X ( REAL ( IDELE(NN) ) ) ) 

IRJ=( IFIX( AIMAG( IDELE(NN) ) ) ) 
CUR1=CABS(FCUR1 ) 

CUR2=CABS(FCUR2) 

IF(CUR1 .LT. 0.001 ) GO TO 400 
WR I TE ( 4 1 , 350)1 FLT, NN, CUR1 , CUR2 
3 50 FORMAT( 2 1 , 2F ) 

400 CONTINUE 

500 CONTINUE 

DO 600 1= 1 , IBB 
READ( 4 2 , 525) IL, JL, ALCUR 
IF( IL. EQ. 0)GO TO 610 
IF(REAL( ALCUR) .LE. 0 .0)GO TO 600 
BLCUR=CABS( ALCUR) 

CLCUR=FIPU*BLCUR*RKVBA/EIB(JL) 

BLCUR = FI PU*BLCUR*R1CVBA/EIB( IL) 

WRITE(43 , 525) IL, JL, BLCUR, CLCUR 
FORMAT ( 21 , 2F) 



5 25 
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6 00 
610 

60 2 



CONTINUE 
WKITE( A 1,602 ) 

WRITE( A3 , 602 ) 

FORMAT ( ' 0 0 0 0 0 ' ) 

CLOSE(UNIT=36) 

CLOSE ( UN1T=37 ) 
CLOSE(UNIT=38) 

CLOSE ( UNIT=39 ) 
CL0SE(UNIT=42) 
CL0SE(UNIT=43) 

STOP 

END 
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**************************************************** 

C TITLE : COORD. FA ; USES DATA FROM 1 NPUT . F A , F ALT . F A TO 

COORDINATE' 

C THE PROTECTION DEVICES FOR THE WHOLE DISTRIBUTION 

SYSTEM. 

*************************************************** 



> 



OPEN( UNIT=30, FILE= ' RE LA. DAT' ) 

OPEN(UNIT=31 ,FILE=' LOCUR.DAT' ) 

OPEN( UNIT=3 2, F I LE=' LCUR.DAT' ) 

OPEN( UNIT=3 3 , FILE= ' FCUR. DAT' ) 

OPEN( UNIT=34 , FILE= ' TEST. DAT' ) 

0 P E N ( U N I T = 3 5 , FILE= ' XFOR. DAT' ) 

OPEN( UN I T = 3 6 , FILE= ' DEV .DAT' ) 

OPE N( UNI T=3 7 , FI LE= ' ADEV.DAT' ) 

DIMENSION EB( 30) , CUREL( 200) , CURLO( 30) , Z( 100) 

1 , AIDEVL(30,4) , A I DE V ( 100,7 ) , XF R( 2 00 , 2) , IID( 3) , RID( 100 

2 , DEV( 20 , 0/4 ) 

COMPLEX Z , AIDEV, Cl , C2 , Z1 , Z2 , RID, CUREL, CU1 , CURLO, XFR 



**************************************** 



******** 



* 



C READ IN LOAD AND BUS VOLTAGE DATA 

*************************************************** 

READ( 30 , 5 ) R 
5 FORMAT(F) 

READ( 34 , 1 0 ) I B , E 1 

10 FORMAT(I,F) 

11 READ( 30 , 12)1 

12 FORMAT(I) 

IF(I.NE.O)GO TO 11 

1 3 READ( 30 , 14 ) I , El 

14 FORM AT ( I , F ) 

IF( I . EQ .0) GO TO 15 
EB( I ) = E 1 
GO TO 13 

15 DO 20 1 = 1 , IB 
READ( 3 1 , 1 6 ) J , CU2 

16 FORMAT( I , F ) 

IF(J.EQ.O)GO TO 21 

CURLO( J)=CMPLX( CU2 * 1 . 25, 0.0) 

20 CONTINUE 

21 DO 24 1=1 , 20 
READ(36,22)K,(DEV(I,J),J=0,4) 

IF(K.EQ.O)GO TO 23 

2 2 FORMAT ( I , 5 F ) 

24 CONTINUE 

************************************************* 



C NOW TO SET ALL LOAD DEVICES 

*************************************************** 

23 TYPE 25 



121 



25 FORMAT( ' ALL LOAD CURRENTS ARE READY, USE DEVICE 

NUMBER LIST,') 

DO 100 1=1 , I B 

I F ( CABS(CURLO( I) ) .EQ.O)GO TO 100 
1 1)= I 

30 TYPE 35, ID 

35 FORMAT( ' SELECT THE DEVICE FOR THE LOAD AT 

BUS' , 13, ' . . , '$) 

ACCEPT 40 , I DEV 
AIDEVL(ID, 1 ) = FLOAT( I DEV) 

40 FORM AT( I ) 

CURLI=REAL( CURLO( I ) ) 

CALL 

SELECT( CURLI , EB( I) , AIDEVL( I , 1 ) , AIDEVL( I , 2) , AIDEVL( I , 3) 

1 , AIDEVL( 1,4)) 

100 CONTINUE 

********************************************** 

C NOW TO BEGIN COORDINATING ELEMENTS IN THE 

DISTRIBUTION 

*********** **************** ************** ******* 

110 REA I) (34, 12Q)NN,I,J, K , R , X 

120 FORMAT( 4 I , 2F) 

I F ( N N . E Q . 0 ) G 0 TO 130 
NE= NN 
CI = I 
CJ = J 

R I D ( N N ) = C M P L X ( Cl , C.l ) 

Z( NN ) = CMP LX ( R , X) 

GO TO 110 

130 DO 140 I E = 1 , NE 

READ ( 3 2 , 135) 1 , J,CU1 

135 FORMATC 21 , 2F) 

I F ( I . E Q . 0 ) G 0 TO 140 
I A = 0 

CALL SERCH(NE, RID, IE, I , J, I A, I ID) 

Z1=0.0*Z1 
Z2 = 0 . 0*Z2 

DO 136 I E L= 1 , I A 
Z 1 = 1 /( Z2 + Z( I I D( I EL) ) ) 

Z 2= 1 / Z 1 

136 CONTINUE 

DO 138 I EL= 1 , I A 

CUREL( I I D ( IA) ) = (Z2)*( 1 / Z ( IID( I EL) ) ) * C U 1 
138 CONTINUE 

IF( IA.GT. 1)TYPE 101, J 

101 FORMATC' CAUTION THERE IS MORE THAN ONE 

SOURCE 

1 0 F POWER TO THE BUS' 14) 

140 CONTINUE 
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1 A 5 



147 



DO 150 



I A= 0 



I = 1 , N E 

READ( 35 , 14 5 ) I J, J , TINKSH, TUST1 , TWST2,TS1 
FORMAT( 21 , 4 F ) 

I F ( I J . EQ . 0) GO TO 151 

CALL SERCH( NE, RID, I , I J , J , I A, 11D) 

I A = 1 

XFR( IID( I A) , 1 )= CMPLX (TS1 , TINRSH) 
XFR(IID(IA) , 2)=CMPLX(TWST1 , TWST2) 



150 CONTINUE 

************** ****************** ************* 



C NOW TO DETERMINE WHAT ELEMENTS HAVE LOAD BUSSES 

CONNECTED 

**************************************************** 

151 DO 190 1=1 , NE 

1 R= 1 F I X ( RE AL( R I D( 1 ) ) ) 

IM=IFIX( AIMAG( R1D(1))) 

IF( AIDEVL( IR, 1) .NE.O) GO TO 155 
IF( AIDEVL( IM, 1 ) • NE.O)GO TO 156 
GO TO 190 

155 11= IFIX( REAL( RID( I ) ) ) 

JJ=IF1X( AIMAG( RID( I ) ) ) 

GO TO 157 

156 11 = IFIX(AIM AG (RID(l))) 

JJ=IFIX( REAL( RID( I ) ) ) 

157 CII=II 

CJ J=J J 

A I D E V ( I , 1 ) = CMPLX( CJJ , Cl I) 

A1 DEV ( I , 2 ) = CMPLX ( 0 . 0 , AI DEVL( 11,1)) 

AI DE V ( I , 4 )= CMPLX ( AIDE VL( II,2)»AIDEVL(II,3)) 
IF( AIDEVL(1 1 , 1 ) . GE. 35) A I DEV ( I , 5)=CMPLX(0 .0 , 

1 AI1)EVL( 11,4)) 

IF(A1DEVL(1I,1).LT.0)A1DEV(1,7)=CMPLX(0.0,AIDEVL(II > 4 

)) 

159 TYPE 158, JJ, II 

158 FORMAT ( ' SELECT DEVICE FOR 

BUS' , 13 ,' PROTECTING ELEMENT 

1FEEDING'//' LOAD AT BUS ', 14 ,'..,'$ ) 

ACCEPT 160, IDEV 

160 FORMAT ( I ) 

A IDEV 1=FL0AT( IDEV) 

DURR=AIMAG( A I DEV ( 1,4)) 

CURR=REAL( AIDEV( I, 4) ) 

IF( AIMAG( A I DEV ( 1 , 2) ) . GT. 0)CURR=1 0**( DURR) 
CURR=REAL( CUREL( I ) ) 

CALL SELECT( CURR, EB( JJ) , A I DEVI , AIDEV2 , AT. DEV 3, CT1 ) 

CALL 

SETDEV( CUREL, I , J J , II , A I DEV 1 , AIDEV2, AIDEV3 , 

1AIDEV , XFR( I , 1 ) , XFR( I , 2) , RE, EB( J J ) , EB( II) , RID, CT1 , DEV 
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) 

190 

200 

205 
2 1 0 

ELEMENT 

230 



EV ) 

250 

73 , 0 . 0 ) 

EV) 

300 

A 00 
410 



CT1 =0 . 0 

IF( RE. NE. 0)GO TO 159 

CO HTINUE 
I E L = 0 

NEL = NE- I EL 

IF(REAL(AIDEV(NEL, 1)) .NE.O)GO TO 300 
I N D E X = 1 

I=IFIX(I!EAL(RID(NEL) ) ) 

J=IFIX(AIMAG(RID(NEL) ) ) 

IF(AIMAG(AIDEV(NEL,2) ) . NE . 0 ) GO TO 250 
TYPE 210,J,NEL 
RE = 0 

FORM AT ( ' ENTER THE DEVICE FOR THE BUS', 13,' OF 
,14) 

ACCEPT 230,1 DEV 
FORMAT( I ) 

IF(IDEV.EQ. 100)GO TO 400 
AIDEV 1 =FLOAT( 1DEV) 

CURR=AIMAG( CUREL( NEL) ) 

CALL S ELECT ( CURR, EB( j) , AIDEV 1 , AIDEV 2 , AIDEV3,CT1) 

IF(AI l)E VI. GT.34)AIDEV(NEL,5) = CMPLX(0.0, .414973) 

CALL SETDEV( CUREL, NEL, I , J , AIDEV 1 , AIDE V2 , AIDEV 3, 

1 AIDEV , XFR( NEL, 1 ) ,XFR( NEL, 2) , RE, EB( I ) ,EB( J) , RID, CT1 , D 

I F ( RE.NE.O)GO TO 205 
TYPE 210, I, NEL 
ACCEPT 230 , 1 DEV 
AI DEV 1 = FLOAT ( IDEV) 

RE = 0 

CURR= REAL (CUREL (NEL) ) 

CALL SELECT( CURR, EB(I), AIDE VI, AIDEV 2, AIDEV3.CT1) 

IF( AIDEV1 . GT . 34 ) AIDEV( NEL, 5 )=AIDEV( NEL, 5 )+CMPLX( . 4 149 

CALL SETDEV( CUREL, NEL, I, J, A I DEVI , AI DEV 2 , AI DEV 3 , 

1 AIDEV, XFR( NEL, 1),XFR(NEL,2),RE,ER(1),EB(J),RID,CT1,D 

IF(RE.NE.Q)GO TO 250 
I E L= I E L+ 1 

IF( IEL. LT. NE) GO TO 200 

IF( INDEX. EQ.O)GO TO 400 

IE L=0 

INDEX= 0 

GO TO 200 

TYPE 410 

PRINT 410 

FORMAT( ' DEVICE LISTING') 

DO 440 1=1 , NE 

TYPE 420,1, (AIDEV( 1,J) ,J=1 ,7) 
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PRINT 420,1, (AIDEV(I,J),J=1,7) 

WRITE ( 3 7 , 42 0) I , ( AIDEV( I , J) , J=1 , 7 ) 

4 20 F ORMAT ( 1 , 1 4 F ) 

440 CONTINUE 

W R I T E ( 37 , 442 ) 

442 FORMAT ( ' 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 

0.0 

1 0.0 0.0 0.0' ) 

2000 CLOSE( UNIT=30) 

CL0SE( UNIT=3 1 ) 

CLOSE(UNIT=32) 

CLOSE ( UNIT=33 ) 

CLOSE ( UNIT=34 ) 

CL0SE(UNIT=36) 

CLOSE( UMIT = 3 7 ) 

STOP 

END 

************************************************** 

************************************************** 

C SUBROUTINE S E RCH ; F I N DS THE ELEMENT CORRESPONDING TO 

TWO 

C BUSSES. 

************************************************** 
SUBROUTINE SERCH( NE , RID,NN,I,J,IA,IID) 

DIMENSION RID( 100) ,IID(3) 

COMPLEX R I D , C 1 , C 2 
I A = 0 

DO 30 I E L = 1 , 100 

CJ = J 
CI = I 

C 1 = CMP LX (Cl , CJ) 

C2=CMPLX( CJ , Cl ) 

I F ( RID( IEL) . EQ. Cl )G0 TO 20 
IF(RID(IEL).EQ.C2)G0 TO 20 
GO TO 30 

2 0 I A= I A+ 1 

1 1 D ( I A ) = I E L 

30 CONTINUE 

31 RETURN 
END 

*************************** * ****** **************** 
************************************************ 

C SUBROUTINE REFALT; READS FCUR.DAT TO OUTPUT FAULT 

CURRENTS 

C AT THE BUSSES 

************************************************* 

SUBROUTINE REFALT ( I FALT , NEL , F ALT , CU RE L ) 

10 R E A D ( 3 3 , 15)1 EF, NN.TFAL, SFAL 

15 FORMAT (21, 2 F ) 
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I F ( IEF.EQ.O)GO TO 20 
IF ( I EF. ME. IF ALT) GO TO 10 
IF( NN. NE . NEL) GO TO 10 
F ALT 1 = TF AL 
F ALT2 = S F AL 

F ALT= AMI N 1 ( F ALT 1 , FALT2) 

FL0CUR=3 . 0*CUREL 

IF( FALT.LE . F LOC UR ) F ALT= AMAX 1 ( F ALT 1 , FALT2) 
FALT= AMAX 1 ( F ALT 1 , F ALT 2 ) 

20 REWIND 33 

RETURN 



END 

i: k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 
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C SUBROUTINE SETDEV; SETS THE DEVICE INSTANTANEOUS 

ELEMENT 

C DEPENDING ON THE FAULT CURRENT AT THAT BUS 

DOES 

C COORDINATE WITH OTHER DEVICES 



* k k * k k k * k k * k 
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SUBROUTINE 

SETDEV ( CUR EL, I EL, I , J , AIDEV 1 , AIDEV2 , AIDE V 3 , AIDEV 
1 , XFR1 , XFR2 , RE, EB1 , EB2 , RID, CT, DEV) 

DIMENSION 

A I DEV ( 10 0,7) , CUREL( 100) , DEV (20, 0/4) , RID( 100) 

COMPLEX AIDEV, XFR1 , XFR2, CUREL, RID 
IDEV1=IFIX(AIDEV 1 ) 

30 IF( IDEV1 .LT.0)CALL 

SETREL(CUREL, IEL, I, J, AIDEV, AIDEV1 , 

1AIDEV2 , A I DEV 3 , XFR1 , XFR2, RE, EB1 ,EB2 , RID, CT, DEV) 

IF (IDE VI .LT.O)GO TO 100 
IF( IDEV1 .EQ. 99 ) GO TO 100 

IF( I DEVI .LT. 3 5 . AND. I DEVI . GT . 0 ) CALL SETFUS ( CUR EL ( I EL 
1 ) , I , J, IEL, AIDEV, I DEVI , AIDEV 2 , AIDEV 3 , XFR1 , XFR2 , RE, EB1 

,EB2, 

2 DEV, RID) 

IF( I DEVI .GE. 35. AND. I DEVI .LT. 99) CALL 
SETCKB( CUREL( IEL) , 

1 I , J , IEL, AIDEV, I DEVI , AIDEV2, AIDEV 3 , CT, EB1 , EB2) 

100 RETURN 

END 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



C SUBROUTINE S E TF U S ; C OORD I NAT E S AND SELECTS A FUSE 

****** A * * * * * *********** * * * * * * ******************* * 



SUBROUTINE 

SETFUS( CUREL, I , J , IEL, AIDEV , I DEV 1 , A I DEV 2 , A I DEV 3 
1 , X F R 1 , X F R 2 , R E , E B 1 , E B 2 , DE V , R I D ) 
DIMENSION 

A I DEV ( 100 , 6) , CUREL( 100) , DEV( 20, 0/A ) , RIP( 1 00) 



126 



COMPLEX A I DEV , XFK1 , XFR2 , CUREL, RID 
JJ=IFIX( AIMAG( AIDEV( I EL, 2) ) ) 

I F ( J J . EQ. 0) CALL XFDEV ( 1 , J , I E L , A I DEV , COOF ) 

I F( J J . NE . 0) CALL XF DEV ( 2 , J , I EL , AI DE V , COOF ) 

IF( JJ. EQ. 0) CURR= AIM AG ( CIJKEL ( I EL) ) 

1 F ( J J . NE. 0) CURR=REAL( CUREL ( I EL) ) 

IF( J J . EQ . 0 . AND. COOF . LT . CURR) COOF=CURR 

IF( JJ.NE.O)COOF=(EB2/EB1 )*COOF 

I F ( JJ . NE .0 . AND. COOF . LT . CURR) COOF= CURR 

IF( IDEV1 .LE. 15) CALL SFUS E ( COOF , I DEV 1 , OF SET , OF SET 1 ) 

IF( IDEV1 .CT. 15) CALL S BO R FU ( CO OF , I DE V 1 , OF SET , OF SET 1 ) 

IF(JJ.EQ.O) GO TO 20 

CI = I 

CJ = J 

AIDEV( I EL, 1 ) = CMPLX( Cl , CJ ) 

A I DEV ( I EL, 2)=AIDEV( I EL, 2) + CMPLX( FLOAT( I DEVI ) ,0.0) 

A I DEV (I EL, 3) = CMPLX( OF SET, 0FSET1 ) 

GO TO 25 
20 CI=I 

CJ = J 

A I DEV ( I EL, 1 ) = CMP LX ( C I , CJ) 

AI DEV ( I EL, 2 ) = CMP LX ( 0 .0,FL0AT(IDEV1)) 

AIDEV( I EL, 4)=CMPLX( OF SET, OF SET 1 ) 

GO TO 50 

25 IF( CABS( XFR1 ) .EQ.O)GO TO 50 

AX= ALOG 1 0 ( (AIMAG(XFRl))/3.0) 

IF(REAL(AIDEV(IEL, 2) ) .LE. 15)G0 TO 35 
IF( OFSET.LE. AX ) CU R= 1 0 * * AX 

IF( OF SET. LE. AX) CALL SBORFU( CUR, IDEV,0FSET,0FSET1) 

GO TO 40 

35 IF (OFSET.LE. AX ) C U R= 1 0 * * AX 

IF( OFSET.LE. AX) CALL SFUS E (CUR, I DEV, OF SET, OF SET 1 ) 

4 0 AIDEV( I EL, 2 ) = CMPLX( FLOAT( I DE V ) , AI M AG ( AI DE V ( I EL , 2 ) ) ) 

A I DEV ( I EL, 3)=CMPLX( OF SET, 0FSET1 ) 

AR=ALOG10(REAL(XFR2) ) 

CALL FU SE ( 2 , REAL ( AI DE V ( IEL,2)),AR,0FSET1,TVAL1) 
TVAL=1 0**TVAL1 
IF(TVAL.GT.2.0)GO TO 43 
GO TO 50 
43 TYPE 45 

45 FORMAT( ' UNABLE TO PUT FUSE IN THIS LOCATION 

SUBSTITUE 

1 ANOTHER DEVICE') 

RE= 1 .0 
GO TO 300 

50 CALL REFALT(J, IEL, FALT, CUREL( I EL) ) 

PAUSE 55 
I F ( J J . EQ. 0) CALL 

5 F ALT 1 ( A I DEV , CURE L , R I D , I E L , I , J , T V AL 1 , DE V , 
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1CVAL) 

IF( J J . NE . 0)CALL S F ALT 2 ( A I DE V , I EL , I , J , TVAL 1 , DEV , 

1 F ALT , CVAL) 

PAUSE 60 

ADEV= FLOAT ( IDEV1 ) 

55 CALL FUSE( 1 , ADEV, CVAL, OFSET, TVAL) 

IF( TVAL. GT. TVAL1 )GO TO 200 
ADEV= ADE V+ 1 
NEGDEV=-IFIX( ADEV) 

I F ( ADEV.LE. 1 5 ) C ALL S FU S E ( C OOF , NEG DEV , 0 F SET , OF SET 1 ) 

1 F ( ADEV.GT. 15)CALL SBORFU ( C OOF , N EGDE V , 0 F S ET , OF S ET 1 ) 

GO TO 55 

200 IF( JJ.NE.O)GO TO 250 

AI DEV ( I EL , 2 ) = CMP LX ( 0 . 0 , ADEV ) 

A I DEV ( I EL, 4 ) = CMP LX (OFSET, OFSET 1 ) 

GO TO 300 

250 A I DEV ( I EL, 2) = CMPLX( ADEV , AIM AG ( A I DEV ( I EL , 2 ) ) ) 

AI DEV ( I EL , 3 ) = CMPLX( OFSET, OFSET! ) 

300 RETURN 

END 

-k v'e & * & * * ******* * * * * *£***/’<&* * *&*£*&£*& >*c ***** * & 
*&****&******&*&>’<***:!:**** ftyc*****^****** ******* *A 

C SUBROUTINE F U S E ; COM P UTE S TIME VALUE FOR A GIVEN LOAD 

CURRENT 

*&&*******&&*&*&*&*&** **&&£&*&&£&£* *&*&&&*****&& 

SUBROUTINE FU SE ( I D , ADE V , CVAL , OF SET , TVAL) 

TV A L=0 . 0 
IDEV=IFIX( ADEV) 

IF( IDEV.GT. 15)GO TO 50 
X= CVAL- OF SET 
I F ( ID.NE. 1 ) GO TO AO 

TVAL=4 .0A80+(-36 . A 2 3 A * X ) + ( 2 2 3 . 7 5 2 6* ( X* * 2 ) ) + ( - 7 5 1 . 279* 

( X* * 3 ) ) 

l+( 1347. 01*(X**4))+(-1222.67*(X**5))+(440. 813 *(X**6)) 
TVAL=TVAL- 2 . 0 
GO TO 100 

AO X= C V AL— OF SET+ . 1 

TV A L=6 . A 2 G+( - 3 7 . 9 2 2*> 1 2 1 . 9 0 9 * ( X** 2 ) ) + (- 2 1 3 . A 20* ( X* 

* 3 ) ) + 

1(201 . 77!*(X**A)) + (- 9 6. 9836 *(X**5)) + (13.60A6*(X**6)) 
TV A L= TV AL- 2 . 0 
I F ( X . G T . 1 . A)TVAL = -1 .95 
GO TO 100 

50 X= CVAL- OF SET 

I F ( ID. NE. 1 ) GO TO 70 

TV A L= 3 . 8 7 6 9 8 + ( - A . 2 6 2 7 6 0 * X ) + ( 1 . 2 1 9 3 A 0 * ( X * * 2 ) ) 
TVAL=TVAL- 2 . 0 
GO TO 100 

70 X=CVAL-OFSET+. 1 
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TVAL-A .61A8A0+(-6.13322A*X)+(3.23181*(X**2))+(-.60571 
36*(X**3) ) 

TVAL=TVAL-2 . 0 
100 RETURN 

END 

& * * * £ * * * * * * * & * * & * £ & * & * * & * * £ £ ft * * * * * * * * * * * * * * 'A * * * & 
*********************************** ****** ******* 

C SUBROUTINE RELAY ; COMPUTE S TIME VALUE FOR A GIVEN 

CURRENT 

C VALUE FOR A SPECIFIC RELAY 

********************************************* 

SUBROUTINE RELAY ( NDEV , CURR, TVAL, TDSET, DEV , AMST) 
DIMENSION DEV(20,0/A) 

TVAL = 0 . 0 
X=AL0C10( CURR) 

IF(X.GT. 1 . 2)X=1 .2 
DO 10 1=0, A 

TVAL=TVAL+DEV(NDEV, I)*(X**I) 

10 CONTINUE 

I F ( NDEV . EQ. 7 ) TVAL=TVAL-2 . 0 
TV AL=TVAL+ TDSET 

I F ( CURR. GT. AMST . AND . AMST. ME. 0)TVAL = -2 . 0 

RETURN 

END 

************************ * * * * * * * * * * * * * * * * * * * ***** 

C SUBROUTINE CKTBKR ; COMPUTES TIME VALUE FOR A GIVEN 

CURRENT 

C VALUE 

***** * * * * * * * * * * ************ * * * ****** * * * * * * * * * * * 

SUBROUTINE CKTBKR ( I D , MDEV , D VAL , OF SET , TV AL , A11VAL) 

CVAL=DVAL-OFSET 

I F ( ID.NE. 1 )G0 TO 50 

CVAL=CVAL+. 1 122 

IF( AMVAL. EQ.0.0)AMVAL=1 0 

IF( CVAL.GT. AM VAL ) GO TO 90 

TVA L= 5 .76908+(-10 . 88261*( CVAL) )+( 23 . 65A 19*( CVAL**2 ) )+ 
( - 3 A . 7 9 2 2 

1*(CVAL**3) ) + ( 33.77 A72*(CVAL**A ) ) + (-21 . 1 A 536* ( CVAL* *5 

) ) + ( 6 . A 

231773*( CVAL**6) ) 

TVAL=TVAL-2 . 0 
GO TO 100 

50 C V A L = C V A L+ . 23 1 5 

AM V A L= A M V A L+ . 0 7 5 
IF( CVAL.GT. AMVAL)GO TO 95 

T V A L = 5 . 85967 2 2+(- A . 39 A 3 3 5 * CVAL ) + ( - . 5 0 6 9 0 2 8 A * ( CVAL* * 2 ) 

) + ( 10 . AO 

12739*(CVAL**3))+(-17.06A956*(CVAL**A ))+(13.0115A3*(C 



VAL**5 ) ) 
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9 0 


2 + (-4 . 010724* ( C V A L * 
TVAL=TVAL- 2 . 0 
GO TO 100 
TVAL = - 2 . 0 


95 


GO TO 100 

TV A L = ALOG 1 0 ( .025) 


100 


RETURN 




END 



>V ************************ ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft ft 
ftftftft ft ftftftftftft ft ft ftftftftftftftftftftft ft ft ftftftftft ft ft ftftftftftftftft ft ftftftft ft 

C SUBROUTINE SETREL ; COORD INATES AND SELECTS A RELAY 

******** * *************************************** 

SUBROUTINE 

SETREL(CUREL,IEL,I,J,AIDEV,AIDEV1,AIDEV2, AIDEV3 , 

1XFR1 , XFR2 , RE, EB1 , EB2 , RID, CT, DEV) 

DIMENSION 

A I DEV ( 100 , 7 ) , IIDD(30) , RID( 100) , IID( 3) , CUREL( 1 00) 

1 , DEV (20, 0/4) , FT( 30 , 2) 

COMPLEX AIDEV, CUREL.XFR1 ,XFR2,RID 

COOF = 0 . 0 

CURR=0 . 0 

CDUR = 0 . 0 

TYPE 5 

5 FORMAT( ' DOES THIS RELAY HAVE AN INST. ELEMENT, IF 

YES 0' ) 

ACCEPT 7 , I L 
7 FORMAT(I) 

JJ=IFIX(AIMAG(AIDEV( IEL, 2) ) ) 

MD KV=IFIX(-AI DEV 1 ) 

IF(JJ.EQ.O)AIDEV(IEL, 1 ) = CMP LX ( F LO AT ( I ) , FLOAT ( J ) ) 
IF(JJ.EQ.0)AIDEV(IEL, 2 ) = CM P LX ( 0 . 0 , AIDEVl ) 

IF( JJ.NE .Q)AIDEV(IEL, 2)=CMPLX( AIDEVl , 0 .0)+AIDEV( IEL, 2 

) 

IF( J J . EQ. 0) CALL XFDEV( 1 , J , IEL, AIDEV , COOF) 

IF( JJ.NE. 0) CALL XFDEV ( 2 , J , I E L , A I DE V , COOF ) 

I F ( JJ. EQ. 0)CURR = AIMAG( CUREL( IEL) ) 

IF( J J . NE .0)CURR = REAL( CUREL( IEL) ) 

IF( JJ . EQ. 0 . AND . COOF . LT. CURR) COOF=CURR 
IF( JJ . NE. 0)COOF=( EB2/EB1 ) *COOF 
CDUR=ALOG 1 0( COOF) 

IF( J J . NE . 0 . AND . COOF . LT. CURR) COOF= CU RR 
IF( CDUR.GT.AIDEV2)AIDEV1 = 4 .0 
IF( CD UR. GT. AIDEV 2 . AND. J J . EQ. 0) CALL 
SELECT( COOF, EB2 , AIDEV 1 , 

1 AIDEV2 , AIDEV3 , CT) 

I F ( CDUR. GT. AIDEV 2 . AND. JJ . NE . 0)CALL 
SELF.CT(COOF,EBl , AIDEVl, 

1 AI D E V 2 , AIDEV3 , CT) 

CURL=AIMAG( CUREL( IEL) ) 
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8 

10 



)))) 



IF( JJ.NE.0)A1DEV(1EL, 3) = CMPLX( A I DEV 2, AIDEV3) 
IF(JJ.EQ.0)A1DEV(1EL,4) = CM1 > LX(AIDEV2,AIDEV3) 

1 F ( J J . NE . 0)CURL=REAL( C U R E L ( 1EL) ) 

IF( J J . NE.Q)A1DEV( 1EL, 7 ) = CMPLX( CT, 0 .0) + AlDEV( I EL, 7 ) 

IF( J J . EQ. 0)AIDEV( IEL, 7 ) = CMPLX( 0 . 0, CT) 

1F( 1L. NE.O)GO TO 100 

CALL REFALT(J, IEL, F ALT 1 , CURL) 

IF(JJ.NE.O ) GO TO 32 
FALT1=FALT1*(EB1 /EB2 ) 

FALT 1 = F ALT 1 / ( 1 0 * * ( R E AL ( A 1 DE V ( I EL, 4) )) ) 

1N = 0 

DO 10 JI=1 , 30 

CALL SERCH( 1 ,R1D, 1EL,J,J1, I A , 1 I D ) 

IF( IA.EQ.O) GO TO 10 
DO 8 J A= 1 , IA 

IN=IN+1 

1 1 1) D ( IN) = 1ID( 1A) 

CONTINUE 

CONTINUE 
DO 30 I N N = 1 , IN 

ID=IIDD( INN) 

IF(ID.EQ.IEL) CO TO 30 
I J=IFIX( AIM AG ( RID( I D) ) ) 

CURJ=AIMAG(CUREL(ID) ) 

CALL REFALT( I J , I D, FALT , CUP.J ) 
FALT2=ALOG10(FALT) 

AID=REAL(AIDEV( ID, 2) ) 

IA1D=IFIX(AID) 

I F ( IAID.GT.O)GO TO 16 

I F ( IAID.LT. 0)FALT0 = FALT/ ( 10**(REAL(AIDEV( ID, 3 



N DE V = - I F I X ( REAL ( A I D E V ( ID, 2) ) ) 

K F A L = I F I X ( 1 0*F ALTO ) 

DO 15 K=1 5 , KFAL, 1 0 
RK=K* . 1 

IF(ALOG10(RK) .GT. 1 . 2) GO TO 16 
CALL 

RELAY ( N DEV , RK, TVAL1 , RE AL ( A I DE V ( I D , 6 ) ) 

1 , DEV , REAL(AIDEV( ID, 5) ) ) 

CALL 

RELAY (MDEV , RK, TVAL2 , TM SET, DEV, 0 .0) 

TIMEO=( 1 0**TVAL1 ) + . 4 

TIMER=TVAL2 

T1MES=ALOG10(TIMEO) 

TMSET1 =TIMES- TIMER 
IF(TMSET1 .LE.O)GO TO 15 
TM S ET=TMS E T 1 + TM S ET 



1 5 

1 6 



CONTINUE 
I F ( IAID.LT. Q)GO 



TO 2 5 
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OF SET = A 1M AG (All)EV(lD, 3 ) ) 
FNSET=FALT2-OFSET 
IF(FNSET.GT. 1 . 5 ) FALT2=OF 5ET+ 1 .5 
KV AL= I F 1 X ( 10*0 F SET) 

KF AL= 1 F I X ( 10* F ALT 2 ) 

DO 20 K = KVAL , KF AL 



RK=K* . 1 

IF( AID. LT. 35) CALL FUS E ( 2 , A 1 D , RK , OF S ET , TVA L 1 ) 

I F ( AID. GE . 35) CALL 

CKTBKR ( 2 , I FIX (A ID) , RK , OF SET, TVAL1 , R E A L ( 

1 A I DEV ( IEL, 5) ) ) 

CUR=( 10**RK) / ( 10**(REAL( A I DEV ( I EL, 3) ) 



)) 



CALL 

RELAY (MDEV , CUR, TVAL2 , TMSET, DEV , 0 . 0 ) 

T I M E 0 = 1 0**TVAL1+. 4 

TIMER=TVAL2 

T IHES= ALOC 1 O(TIMEO) 

TMSET 1 = TIMES- TIMER 
IF(TMSET1 .LE.O)GO TO 20 
TMSET=TMSET+TMSET1 

20 CONTINUE 

25 GO TO 30 

30 CONTINUE 

A I DEV (IEL, 6 ) = CMPLX( 0.0, TMSET) 

TMS ET = 0 . 0 
GO TO 45 

32 AID=AIMAG( AIDEV( IEL, 2) ) 

FALT2=ALOG10( FALT1 ) 

IF(AID.GT.O)GO TO 35 
FALT1=FALT1*(EB1/EB2) 

IF( AID. LT. 0)FALT=EALT1 / ( 1 0 * * ( RE AL( A I DEV ( I E L , 4 ) ) ) ) 
NDE V = - I F I X ( AIM AG ( A I DEV (IEL, 2))) 

KFAL=IFIX( FALT* 10) 

DIF = KF AL- 1 5 
ID=IFIX(DIF/50) 

I F ( I D . LE . 0 ) I D= 1 
DO 3 3 K= 1 5 , KF AL , I D 
RK=K* . 1 

RK2 = RK*( EB2/EB1 ) * ( 1 0** ( REAL( A I DEV ( I EL , 4 ) ) ) / 1 0** 

1 (REAL( AIDEV( IEL, 3) ) ) ) 

IF( RK. GT. 20)00 TO 34 
CALL 

RELAY (N DEV , RK, TVAL1 , AI MAG ( AI DE V ( I EL , 6) ) , DEV , 

1AIMAG( AIDEV( IEL, 5) ) ) 

CALL RELAY ( MDEV , RK 2 , T V A L2 , TM S ET , D E V , 0 . 0 ) 

T 1 M E 0 = 1 0**TVAL1 + . 4 
TIMES=AL0G1 O(TIMEO) 

T1MER=TVAL2 



I 
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TM S E T 1 = T I M E S - T I M E R 
IF(TMSET1 .LE.O)GO TO 33 
TM SET= TM SE T 1 + TM S E T 

33 CONTINUE 

3 4 AIDEV ( 1EL, 6)=AIDEV( IEL, 6) + CMI’LX( TMSET, 0 .0) 

35 TMSET=0.0 

IF(AID.LT.O)GO TO 45 
OFSET = AIMAG( AIDEV (I EL, 4) ) 

FNSET=FALT2-OFSET 
IF(FNSET.GT. 1 . 5 ) F ALT2 = OF S ET+ 1 .5 
KVAL=IFIX( 1 0* OF SET) 

KF AL= I F I X ( 1 0* F ALT2 ) 

DO 40 K = K VAL , KF AL 
RK = K* . 1 

I F ( AID. LT. 35 ) CALL FUSE ( 2 , AID , RK , OF SET , TVAL 1 ) 

I F ( AID. GE . 35)CALL 

CKTBKR( 2 , IFIX(AID) , RK , OF SET , TVAL 1 , AI MAG ( 

1 AIDEV ( IEL, 5) ) ) 

CUR = RK-REAL( AIDEV (IEL, 3) ) 

CURTA= 1 0* * CU R 

CUR=( ER2/EB1 )*( 1 0**RK) 

CURTA=CUR/ ( 1 0** ( RE AL( AIDEV ( IEL, 3 ) ) ) ) 
IF(CURTA.LT. 1 . 5 ) GO TO 40 

CALL R ELAY ( MDEV , CURTA, TVAL2 , TMS ET , DEV , 0 . 0 ) 
TIMEO=l 0**TVAL1+. 4 
TIMER= TVAL2 
TIMES=AL0G1 O(TIMEO) 

TMS ET 1 = T IME S- T IME R 
IF( TMSET1 .LE. 0)GO TO 40 
TMSET=TMSET+TMSET 1 

40 CONTINUE 

AIDEV ( IEL, 6)=AIDEV( IEL, 6 ) + CMP LX ( TMS ET , 0 .0) 

45 TMSET = 0 . 0 

IF( J J . EQ. 0 . AND. F ALT 1 . GT. 6)FALT1 = 6 

I F ( IL.EQ.O. AND. JJ.EQ.O ) AIDEV (IEL, 5 ) = CMP LX ( 0 . 0 , F ALT 1 ) 
FALT 1 = FALT 1 / ( 1 0* * REAL ( AIDEV ( IEL, 3) ) ) 

I F ( F ALT 1 . GT . 6 ) F ALT 1=6 

IF( IL. EQ . 0 . AND. JJ . NE. 0) AIDEV ( IEL, 5) = AIDEV( IEL, 5 ) + CMPL 

X ( F ALT 1 , 

10.0) 

100 RETURN 

END 

A*********************************************** 

C SUBROUTINE SF ALT 1 ; PRO V I DE S THE DEVICE WITH THE 

LONGEST 

C TIME AT FAULT, FIRST DEV IN ELEMENT. 

*********************************************** 

SUBROUTINE SF ALT 1 ( A I DEV , CUR EL , R I D , IEL, I , J , T VA L , DE V , 

1C VAL) 
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DIMENSION AIDEV( 100,6) ,CUKEL( 100) , R I D ( 100) , I I D ( 3 ) , 

1 1 1 DD ( 3 0 ) , DEV (20,0/4) 

COMPLEX AIDEV, CUREL, RID 
DO 10 J 1= 1 , 30 

CALL SERCH( 1 , RID, ILL, J , JI , I A, I ID) 

IF( IA.EQ.O)GO TO 10 
DO 8 J A= 1 , I A 

I N= I N+ 1 

IIDD(IN) = IID( IA) 

8 CONTINUE 

10 CONTINUE 

DO 30 INN-1, IN 

ID=IIDD( INN) 

IJ=IFIX(AIMAG(RID( ID) ) ) 

CURJ=AIMAG( CUREL( ID) ) 

CALL REF ALT( I J , ID, FALT, CURJ ) 
FALT2=AL0G10(FALT) 

AID=REAL(AIDEV(ID, 2)) 

I A I D= I F IX(AID) 

IF( IAID.GT. 32)GO TO 100 
FALT2=REAL( AIDEV ( I EL, 4 ) )+. 5 
FALT0=10**(FALT2) 

IF( IAID.LT. 0)FALT=FALT0/ ( 10** ( R E AL ( A I D E V ( I D , 3 

)))) 

NDEV = -IFIX(REAL( AIDEV ( ID, 2) ) ) 

IF( IAID.LT. 0) CALL RELAY ( NDEV , FALT , TVAL1 , REAL 
1 ( A I DEV ( ID, 6) ) ,DEV, REAL( AIDEV( ID, 5) ) ) 

F AL= REAL (AIDEV(IEL,4 ) )+. 5 
IF( FAL. LT. FALT2)FALT2 = F AL 
IF ( IAID.GT. 0)CALL 
FUSE( 2 , AID, FALT 2 , AIM AG ( AIDEV 
1 ( ID, 3) ) , TVAL1 ) 

IF(TVAL.GE.TVAL1 )GO TO 25 
GO TO 30 

25 IF( IAID.LT. 0 ) CVAL=-F ALT2 

CVAL=F AL 
TVAL=TVAL 1 

30 CONTINUE 

100 RETURN 

END 

********************************************** 

********************************************** 

C SUBROUTINE S F ALT2 ; A S IN 1, COMPUTES TIME FOR LOAD 

SIDE OF 

C ELEMENT AT FALT. 

******************************************** 

SUBROUTINE S FALT 2 (AIDEV , IEL , I , J , T V AL 1 , DEV , F ALT , CVAL ) 
DIMENSION AIDEV( 100,6) ,DEV(20,0/A) 

COMPLEX AIDEV 
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F ALT2 = ALOG 10(FALT) 

AI D = AIM AG ( AI DE V ( I EL , 2 ) ) 

FAL=REAL(AIDEV( I EL, 3) )+.5 
FALT=10**(FAL) 

FALT3 = FALT/ ( 1 0**( AIMAG( A I DEV ( I EL, 3) ) ) ) 

IF( AID.LT.O)CALL 

RELAY (N DEV, FALT3 ,TVAL1 , A I M AG ( A I DE V ( I EL 
1,6)), DEV, A I MAG (A IDE V( I EL, 5) ) ) 

IF( FAL. LT. FALT2)FALT2=FAL 

IF(AID.GT.O)CALL FUSE(2,AID,FAL,AIMAG(AIDEV(IEL,4)) 
1 , TVAL1 ) 

CVAL=F ALT2 

RETURN 

END 



3 0 
32 



E V ( 



DEV ( 



D 



40 

50 



SUBROUTINE XF DEV ( NN, I , I EL , AI DE V , COOF ) 

DIMENSION 1 1 D( 3 0 ) ,AIDEV( 1 00, 6) 

COMPLEX AIDEV 
COOF=0 .0 
IX = 0 

IF( NN. HE . 1 )G0 TO 50 
DO 30 I A= 1 , 100 

AI D= REAL ( AIDEV ( IA, 1 ) ) 

IF( IFIX( AID) .NE. I)GO TO 30 

IX=IX+1 

1 1 D ( I X ) = I A 

CONTINUE 
DO 3 2 K = 1 , IX 
CONTINUE 
DO 40 J A= 1 , IX 

OFSET=AIMAG( AIDEV ( I ID( JA) , 3) ) 

IF(REAL( AIDEV( IID( JA) , 2) ) .LT. 0) OFSET=REAL( AID 
1IID(JA) ,3)) 

I F ( REAL ( A I DE V ( 1 1 D ( J A ) , 2 ) ) . GE . 3 5 ) 0 F S ET= RE AL ( AI 
1IID( JA) ,5) ) + 0F S ET- . 0 5 

IF( AIMAC( AIDEV( I EL, 2 ) ) .LT. 0 . AND. REAL( AIDEV( I I 

1 ( J A ) ,2)) . GT . 0 . AND. REAL (AIDEV ( IID(JA),2)).LT.35) 
20FSET=0F S ET+ . 1 5 

IF( COOF . GT. OFSET) GO TO 40 
COOF=OFSET 

CONTINUE 
GO TO 60 

IF( AIMAG( AIDEV (I EL, 2 ) ) . GT. 0 ) C OOF = A I M AG ( A I DEV( I EL, 4 ) ) 
IF( A I MAG ( AIDEV( I EL, 2) ) . LT. 0 )COOF = REAL( A I DEV ( I EL, 4) ) 

I F ( A I M AG ( A I D E V ( IF.L, 2) ) .GE. 35)COOF = AIMAG(AIDEV( I EL, 5) ) 
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-.075 

1+C00F 

I F ( REAL( AIDEV ( I EL , 2 ) ) . LT . 0 . AN D . AIM AG( A1 DEV ( I EL , 2 ) ) . GT 

. 0 

1 . AND. AIM AG ( AIDEV( I EL, 2) ) .LT. 3 5 ) CO OF = COOF+ . 1 5 
60 COOF=COOF+.025 

D00F=1 0**( COOF) 

COOF = DOOF 

RETURN 

END 

ftftftftftftftftftftftftftftftftftftftftftftftft ft ft ftftftftftftftftftftftftftftftftftftft 
ftft ft ftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftft ft ft ftftftftftftft ft ft 

C SUBROUTINE S ETCKB ; CORD I NATE S AND SELECTS A THERMAL 

MAGNETIC 

C BREAKER. 

ft ftftftftftftftftftftftftftftftftftftftftftftftftftftftft ft ft ftftftftftftftftftftftftftftftft 
ftftftftftftftftftftft ftftftftftftftft ft ftftftft ft ftftftft ftft ft ftft ft ftftftftftft ft ft ft ft ft 

SUBROUTINE S ETCKB ( CUREL , I , J J , I EL , A I DEV , I DEV , OF SET 
1 , OF SET 1 , AMAG.EB1 , EB2) 

DIMENSION AI DE V ( 100,7) ,CON(5) 

COMPLEX AIDEV, CUREL 

CON( 1 ) = 2 .6 

C 0 N ( 2 ) = 3 .5 

CON( 3 ) = 5 . 0 

CON( 4 ) = 6 . 0 

CON(5) = 7 .0 

J=IFIX(AIMAG(AIDEV( I EL, 1 ) ) ) 

I F ( J . EQ . 0 ) CALL XFDE V ( 1 , .1 , I EL , A I DEV , COOF ) 

I F ( J.NE.0)CALL XFDEV(2,J,IEL, AIDEV, COOF) 

IF( J . EQ. 0)DOOF=ALOG1 0( COOF)- AIMAG( AIDEV ( I EL, 5) ) 
IF(J.EQ.O)COOF=10**DOOF 

IF( J.NE. 0)DOOF = ALOG1 0( COOF )- AIM AG ( AIDEV (I EL, 5) ) 

IF( J . NE . 0) COOF = 1 0 * * D 0 0 F 
IF( J . EQ. 0)CURR=AIMAG( CUREL) 

IF( J . NE. 0) CURR= REAL ( CUREL) 

I F ( E B 1 . N E . E B 2 ) G 0 TO 4 5 
IF ( COOF .LT. CURR) COOF =C UR R 
GO TO 46 

45 RE2=(EB2/EB1 ) * C 0 0 F 
IF(RE2.GE. CURR) C OOF =RE2 
IF( RE2 . LT. CURR ) COOF=CURR 

46 CALL SBRKR( COOF , I DEV 3 , OF SET , OFSET1 , DM MY ) 

I F ( IDEV.LT.IDEV3)IDEV=IDEV3 

IF( IDEV.GE. IDEV3)IDEV=IDEV+I 
CI = I 
CJ = J J 

AIDEV (I EL, 1 ) = C M P L X ( Cl , CJ ) 

I F ( J . EQ . 0 ) AI DE V ( I EL , 2 ) = CMP LX ( 0 . 0 , FLO AT( I DE V ) ) 
IF(J.NE.Q)AIDEV(IEL, 2) = CMPLX(FLOAT( I DEV) ,0.0) + A I DEV 
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1 ( I E L , 2 ) 

IF(J . EQ. 0)AIDEV( IEL, 4 )=CMPLX(OFSET, 0FSET1 ) 

IF( J.NE.O)AIDEV(IEL, 3 ) = CMPLX ( OF SET , OF SET 1 ) 

1 F ( J . EQ. Q)FLT=4*AIMAG( CUR EL) 

I F ( J . NE. 0)FLT=4*REAL( CUKEL) 

SET=10**(OFSET— . 12) 

DO 50 K=1 , 5 

TRP=CON(K)*SET 
IF(TRP. GT. FLT)GO TO 60 

50 CONTINUE 

60 IF( AMAG. NE . 0 ) AMAG = ALOG 1 0 ( CON ( K ) ) 

IF( J. EQ. 0) A1DEV( IEL, 5)=CMPLX(0 . 0 , AMAG) 

IF( J.NE.O ) AIDE V( IEL, 5 )= CMPLX ( AMAG , 0 . 0 )+AI DE V ( I E L , 5) 

RETURN 

END 

*****£************&***&&*&*:><*******&*******&** 

C SUBROUTINE SELECT ; SELECTS A PARTICULAR DEVICE 

************** 5'c *************************** 5V 5 V ******* 

SUBROUTINE SELECT(CURLO, EB.ADEVLl , ADEVL2 , ADEVL3 , CT) 
IDEVL1=IFIX( ADEVL1 ) 

54 IF( IDEVL1 .EQ. 0)GO TO 99 
IF( IDEVL1 . EQ. 1 )G0 TO 60 
IF( IDEVL1 . EQ. 2)G0 TO 70 
IF( IDEVLl .EQ. 4)G0 TO 80 
IF( IDEVLl .GT. 4)G0 TO 53 
GO TO 58 

53 TYPE 55 

55 FORMAT(' NO SUCH DEVICE OR DEVICE ERROR, RE-ENTER 

1 DEVICE NR 1-4 . . , ' $ ) 

ACCEPT 56, I DE VL 1 

56 FORMAT(I) 

GO TO 54 

58 IF(EB.GT. 600.0) GO TO 53 

CALL SBRKR( CURLO, IDEVL1 ,ADEVL2,ADEVL3 ,ADEVL4) 
ADEVL1=FL0AT(IDEVL1 ) 

CT= ADEVL4 
GO TO 100 

60 CALL SFUSE( CURLO, IDEVL1 , ADEVL2 , ADEVL3 ) 

ADEVL 1 =FLOAT ( IDEVL1 ) 

GO TO 100 

70 IF(EB.LT.600)GO TO 53 

CALL SBORFU( CURLO, IDEVL1 , AD EV L 2 , ADE VL 3 ) 

ADEVL 1=FL0AT( IDEVL1 ) 

GO TO 100 

80 CALL SRELA( CURLO, IDEVLl , NTAP.NCT, PKUP, CT) 

ADEVL 1=IDEVL1 
ADEVL2 = P KU P 
ADEVL3=NCT* 100+NTAP 
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GO TO 100 


99 


ADE V L 1 = 9 9 




ADEVL2=CURLO 




ADEVL3=CURLO 


100 


CONTINUE 




RETURN 




END 



*********************************************** 



C SUBROUTINE SFUSE;USED TO SELECT A CURRENT LIMITING 

FUSE 

C THIS MAY BE EITHER HIGH OR LOW VOLTAGE 

**************************************************** 



CURLO 

CU 

OG 



1 0 

2* ( DEV 
3 9 * ( 



SUBROUTINE S F US E ( CURLO , NDE V , OF S ET , 0 F SET 1 ) 
IF(NDEV.LT.O)GO TO 10 

DEV = — 1 1 . 1 522 l + ( 23 . 865*ALOG10( CURLO) )-( 1 6 . 4 87 * ( ALOG 1 0 ( 

1 )**2 ) ) + ( 2 . 9 288*( AL0G10( CURLO) **3 ) )+( 3 . 1 9 6 7 2 * ( ALOG 1 0 ( 

2RLO)**4))-(l .48877*(AL0G10(CURL0)**5) ) + ( . 1776431*(AL 

3 1 0 ( CURLO) * *6 ) ) 

NDEV=IFIX(DEV)+1 

IF( MDEV.GT. 1 5)NDEV=1 5 

DEV=ABS( FLOAT(NDEV) ) 

OFSET=. 953 2- ( . 04 64 3 3*DEV )+ ( . 1 1 7 88 5 * ( DEV** 2 ) )+(- .02935 
l**3))+(.00337*(DEV**4))+(-.0001847*(DEV**5))+( .00000 

2 D E V * * 6 ) ) 

0FSET1 =OFSET+. 1 
IF(NDEV.LT.-15)NDEV=— 15 
IF( NDEV.GT. 1 5)HDEV=1 5 
RETURN 

END 



*************************************************** 



C SUBROUTINE SBORFU;USED TO SELECT A BORIC ACID HIGH 

VOLTAGE 

C FUSE 

*************************************************** 

SUBROUTINE SBOKFU( CURLO , NDE V , OF SET , OF SET 1 ) 
AL0GC=AL0G10( CURLO) 

I F( NDEV . LT . 0 ) GO TO 10 

DEV=88.1643+(-215.2235*(ALOGC))+(188.225*(ALOGC**2))+ 

(-69. 

1 3463*(ALOGC**3)) + (7.4695*(ALOCC**4 ) )+( 1 . 6 52 5 7 * ( ALOGC 

** 5 ) ) + 

2 ( - . 347786*(ALOGC**6) ) 

NDEV=IFIX(DEV)+1 
I F ( NDEV.GT. 19)NDEV=19. 
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10 DEV=ABS( FL0AT( NDEV) ) 

IF(NDEV.LT.0)DEV=DEV-15 

OFSET=l . 366 1 1 1+. 1 1 34579*DEV+(- . 001 2654*(DEV**2 ) )+( .00 

00213* 

1 (DEV* *3 ) ) 

0FSET1=0FSET+. 1 

IF(NDEV.LT.-19)NDEV=-19 

IF(NDEV.GT. 19)NDEV=19 

I F ( NDEV . LT . 0 ) NDEV=-NDEV+ 1 5 

IF( NDEV . GT.0)NDEV=NDEV+1 5 

RETURN 

END 

**************** ******************************* 

C SUBROUTINE S B RKR ; US ED TO SELECT A THERMAL MAGNETIC 

BREAKER 

*************************************************** 

SUBROUTINE SBRKR( CURLO, NDEV , OF SET , OF SET 1 , AM AG) 
DIMENSION C 0 N ( 5 ) 

CON( 1 )=2 . 6 
CON( 2 ) = 3 .5 
CON( 3 ) = 4 . 6 
C0N( 4 ) = 5 . 6 
C 0 N ( 5 ) = 7 . 0 
TYPE 10 

10 F0RMAT( ' ENTER MANUFACTURER ( N . A . AT PRESENT USE 

FOR FUTURE 

1EXPANSI0N) ' ) 

TYPE 25 

25 FORMAT( ' INPUT TEMPERATURE FOR AMBIENT SHIFT, 1 FOR 

LESS THAN 

110 DEG C , 2 FOR 60 DEG C ,'//'< CR> F OR NO SHIFT..') 
ACCEPT 28 , I OFF 

28 FORMAT(I) 

IF( I OF . EQ . 1 ) OF S = - . 1 
IF(I0F.EQ.2)0FS=.05 
IF(MAN.NE.O)GO TO 80 
IF( CURLO. LE. 600 . )G0 TO 30 
IF( CURLO. LE . 800 .) GO TO 50 
GO TO 70 

3 0 I F ( CU RLO . LT . 3 0 0 . ) MDE V= 1 

IF( CURLO.LT. 300. )0FSET=2. 477 12 

IF ( CURLO .GT. 300 . 0 . AND. CURLO. LE. 350. ) MDE V=2 

IF( CURLO . GT . 300 . 0. AND. CURLO . LE. 350 .) OF SET=2 . 54407 

I F( CURLO. GT. 3 5 0.0 . AND. CURLO . LE . 400 . )MDEV = 3 

I F ( CURLO. GT. 350.0. AND. CURLO. LE. 4 50. )0FSET=2 .60206 

IF( CURLO. GT. 400 .0. AND. CURLO. LE . 450. )MDEV=4 

IF( CURLO. GT. 4 00. 0 . AND. CURLO. LE. 4 50. )0FSET = 2 . 653 2 1 

I F (CURLO. GT. 450.0. AND. CURLO. LF.. 500. )MDEV=5 

IF( CURLO. GT. 450.0 .AND. CURLO. LE. 500. )0FSET=2 .69897 



l 
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IF( CUKLO.GT. 500. )MDEV=6 
I F( CUKLO.GT. 500. ) OF SET=2 . 7 7 8 1 5 
GO TO 30 

50 1F( CURLO. LE. 700 . )MDEV=7 

1F( CURLO. LE. 700 . )0FSET=2 . 8451 
IF( CUKLO.GT . 700. )MDEV=3 
I F ( CUKLO.GT. 700 . )OFSET=2 .90309 
GO TO 80 

70 1F( CURLO. LE. 1000. )MDEV=9 

IF( CURLO. LE. 1 000 . ) OFSET=3 . 0 
I F ( CUKLO.GT. 1000. )MDEV=10 
1F(CURL0.GT. 1000.)0FSET=3.3 
80 HDEV=MDEV+34 

OFSET=OFSET+OFS+. 1122 
OF5ETl=OFSET+. 12 
TYPE 90 

90 FORMAT( ' DOES THE DEVICE HAVE A MAGNETIC 

ELEMENT? (O)YES. ' ) 

ACCEPT 9 5, NOEL 
95 FORMAT(l) 

FCUR=4*CURL0 
S E T = 1 0**(0FSET-. 12) 

DO 100 1=1,5 

TD=SET*CON( I ) 

IF(TD.GT. FCUR)GO TO 105 

100 CONTINUE 

105 AMAG=AL0G10( C0N( I ) ) 

1F( NOEL. NE. 0 ) AM AG=0 . 0 

RETURN 

END 



******* A: *&**&**A::S:A: A* ***&&*£** A: 



******************** 



C SUBROUTINE SRELA;USED TO SELECT A RELAY 



A****************** ********* ********************* 



SUBROUTINE SRELA( CURLO , NDEV , NTAP , NCT, PKUP, CT) 

TYPE 10 

10 FORMAT( SELECT MANUFACTUER (FUTURE EXPANSION)') 

TYPE 20 

20 FORMAT( ' SELECT THE WESTINGHOUSE CO RELAY , REFER TO 

TABLE' , $) 

ACCEPT 25 , NUM 
25 FORMAT(I) 

N D E V = - N U M 

1 F ( CURLO. LT. 30)CT=1 0 

I F( CURLO. GE . 30. AND. CURLO.LT. 1 20) CT= 30 

1F(CURL0.GE. 120)CT=120 

IF(CT.E0. 10) NCT= 1 

I F ( CT. EQ. 30) NCT=2 

1F( CT. EQ. 1 2 0 ) N C T = 3 

TAP=CURL0/CT 



7 1 

72 

73 

74 

75 

76 

77 

78 

79 

8 0 

81 

8 2 



140 



I F ( T A P . LE. 1 . 0 ) G 0 TO 71 

I F ( TAP . LE. 1 . 2) GO TO 72 

I F ( TAP. LE. 1 . 5 ) GO TO 73 

IF( TAP . LE . 2 . 0 ) GO TO 74 

IF(TAP. LE. 2 . 5) GO TO 75 

I F ( TAP. LE. 3 . 0 )G0 TO 76 

I F ( TAP. LE. 3 . 5)G0 TO 77 

I F ( TAP . LE . 4 . 0)GO TO 78 

I F ( TAP. LE. 5 .0)GO TO 79 

IF(TAP.LE.6.0)GO TO 80 

1 F ( TAP. LE. 7 . 0)CO TO 81 

IF(TAP.LE.8.0)GO TO 82 

IF( TAP . LE . 1 0 . 0 ) GO TO 83 

GO TO 84 

NTAP= 1 

TAP= 1 . 0 

GO TO 90 

N T A P = 2 

TAP= 1 . 2 

GO TO 90 

NTAP=3 

TAP= 1 . 5 

GO TO 90 

N T A P = 4 

TAP = 2 . 0 

GO TO 90 

N T A P = 5 

TAP = 2 . 5 

GO TO 90 

K T A P = 6 

TAP = 3 . 0 

GO TO 90 

N T A P = 7 

TAP = 3 . 5 

GO TO 90 

NTAP=8 

TAP=4 . 0 

GO TO 90 

NTAP=9 

TAP = 5 . 0 

GO TO 90 

N T A P = 1 0 

TAP = 6 . 0 

GO TO 90 

N T A P = 1 1 

TAP= 7 .0 

GO TO 90 

N T A P = 1 2 

T A P = 8 . 0 
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GO TO 90 


83 


NTA P= 1 3 . 




T A P = 1 0 . 0 




GO TO 90 


84 


NTAP = 1 4 




TAP= 12.0 


90 


PKU P=ALOG10(TAP*CT) 




RETURN 




END 
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C 

****** 



3 

5 

1 0 

1 5 
20 



25 

30 

40 

45 

50 

55 

60 



70 

75 

80 

400 

410 



420 

440 

445 

450 



TITLE : PLOTD. F4 ; PLOTS DEVICES FROM COORD. F4 
************************************************ 

OPEN(UNIT=30, FILE=' ADEV.DAT' ) 

OPEN( UHIT = 3 1 , FILE = ' DEV . DAT' ) 

OPEH( UNIT=32 , FILE= ' RE LA. DAT' ) 

OPEN( UNIT = 3 3, FILE= ' FCUR. DAT' ) 
OPEN(UNIT=34,FILE='TEST.DAT') 

EXTERNAL T25HAN 

DIMENSION EB( 30) , DEV( 20, 0/4 ) , A I DEV ( 100, 7 ) , RID( 100) 
COMPLEX AI DEV, RID 
READ( 32, 3 ) R 
FORMAT( F) 

READ( 32 , 10)1 
FORMAT( I ) 

IF(I.NE.O)GO TO 5 
READ( 3 2 , 20) I , El 
FORMAT( I , F) 

I F( I . EQ . 0 ) GO TO 25 
EB( I ) =E 1 
CO TO 15 
DO 40 K= 1 , 100 

READ( 30 , 30) I , (AIDEV(I,J),J=1,7) 

FORMAT( I , 1 4 F ) 

I F ( I . EQ . 0 ) GO TO 45 
CONTINUE 

REA D ( 3 4 , 50 ) I B , E 1 
FORMAT( I , F) 

READ ( 34 , 60)NN,I,J,K,R,X 
FORMAT( 41 , 2 F ) 

IF( NN.EQ.O) GO TO 70 
N E = N N 
CI = I 
CJ = J 

RID( NN)=CMPLX(CI , C J ) 

GO TO 55 
DO 80 1=1,20 

READ( 3 1 , 7 5)K, ( DEV( I , J) , J = 0 , 4) 

FORMAT( I , 5 F ) 

IF(K.EQ.O)GO TO 400 

CONTINUE 

TYPE 410 

FORMAT( ' DEVICE LISTING') 

DO 440 1=1, HE 

TYPE 420 , I , ( AIDEV( I , J) , J=l , 7 ) 

PRINT 420,I,(AIDEV(I,J),J=1,7) 

FORMAT ( I , 1 4 F ) 

CONTINUE 
TYPE 450 

F0RMAT( DEVICES HAVE BEEN SET IF YOU WANT TO PLOT 
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GRAPH PAPER 

] ENTER(O) , IF NOT TYPE(l)') 

ACCEPT 456, IP 
456 FORHAT(I) 

IF( IP.NE.O)GO TO 480 
TYPE 460 

460 FORMAT ( ' THIS WILL PLOT CURVES BASED TO THE LOWEST 
VOLTAGE 

1 IN THE SYSTEM. ' ) 

TYPE 461 

461 FORMAT( ' INPUT THE CURRENT OFFSET FOR THE 
GRAPH, INTEGER. . ' ) 

ACCEPT 462 , IAOFF 

462 FORMAT ( I ) 

CALL INITT 
CALL LOG ( IAOFF ) 

CALL FLUSH 

480 TYPE 4 85 

485 FORMAT ( ' INPUT ELEMENT ..') 

ACCEPT 490, IEL 
490 FORMAT ( I ) 

I=IFIX(REAL(RID(IEL) ) ) 

J=IFIX(AIMAG(RID( IEL) ) ) 

EB 1 = EB ( 1 ) 

DO 495 K=1 , IB 

EB 1 = AMI N 1 ( E B 1 , EB(K) ) 

495 CONTINUE 

550 FAC=AL0G10(EB(I)/EB1)-IA0FF 

FAC1=ALOG10(EB(J)/EB1)-IAOFF 
CUR=REAL( AIDEV( IEL, 3) ) 

CUR1=REAL( AIDEV( IEL, 4) ) 

CUR2=AIMAG(AIDEV(IEL, 3) ) 

CUR3 = AIMAG( A I DEV ( IEL, 4 ) ) 

***&********&*****&*&***********************;&******* 

C NOW TO PLOT LOAD SIDE OF THE ELEMENT 

CALL REFALT(J, IEL, FALT, 900.0) 

IF( AIMAG( AIDEV( IEL, 2 ) ) .GT . 0) GO TO 585 
STL=REAL( AIDEV( IEL, 4) ) 

ST P= FALT/ 10** STL 

IF(AIMAG(AIDEV( IEL, 5) ) . NE . 0 ) S TP= AIM AG ( A I DEV ( IEL, 5) ) 
ST1=AL0G1 0( FALT) 

YM=(STL+FAC1) *190+ 175.0 
XM=0 . 0 

CALL MOVE A( XM , YM ) 

STP=10*STP 
ITP=1FIX( STP) 

I F ( ITP.LT. 1 5) TYPE 575 

FORMAT ( ' THE LOAD SIDE RELAY WILL NOT BE VERY 



575 




- 
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SENSITIVE 

1 TO A FAULT' ) 

II)EV = IFIX( A B S ( AIM AG ( AIDEV( I EL, 2 ) ) ) ) 

TYPE 570, I DEV, ITP 
570 FORMAT ( ' I DEV, ' , 21 ) 

D I F= I T P- 1 5 . 

ID=IFIX(DIF/ 200.0) 

IF(ID.LE.0)ID=1 
TVA L3 = 1 0 .0 

TDSET=AIMAG(AIDEV( I EL, 6) ) 

DO 530 K= 1 5 , I T P , I D 

RK= . 1* ( FLOAT ( K ) ) 

CALL RELAY ( I DEV , RK, TVAL, TDSET, DEV ) 
IF(TVAL.CT. TVAL3 )TVAL=TVAL3 
TVAL3=TVAL 

T= 9 00- ( (TVAL+2 . 0)*180) 

R=( ALOG10(RK)+FAC1+STL)* 190+1 75 .0 
CALL DRAWA(T.R) 

580 CONTINUE 

PAUSE 4 

IF(AIMAG(AIDEV( ILL, 5) ) .KE.0)CALL DR AW A ( (900.0) , R) 
GO TO 591 

585 0FSET=CUR1 

M=IFIX(80*OFSET) 

ADEV = AIMAG( AIDEV( I EL, 2) ) 

KIM=IFIX( 100*( ALOGIO(FALT) ) ) 

DO 537 K=M , KIM 

C V A L = K * . 0 1 

IF( ADEV.LT. 35) CALL FUSE(1,ADEV,CVAL,0FSET,TVAL) 

I F ( ADEV . GE . 35)CALL 

CKTBKR( 1 , IFIX( ADEV) , CVAL, OF SET, TVAL, 

1AIMAG( AIDEV( I EL, 5) ) ) 

TVAL=TVAL+2 . 0 

IF(TVAL2.LT. TVAL. AND. ND.NE.O)GO TO 539 

TVAL2=TVAL 

ND= 1 

CALL 

I)RAWA( 9 00- ( ( TVAL) *1 30) , ( CVAL + F AC 1 ) * 1 9 0+ 1 7 5.0) 

537 CONTINUE 

589 OF SET=CUR3 

CALL MOVE A( 0 . 0 , 0 . 0 ) 

N D=0 

M= IFIX(80*0FSET) 

DO 588 K=M , KIM 

C V A L = K * . 0 1 

IF( ADEV .LT. 35)CALL FUSE(2, ADEV, CVAL, OF SET, TVAL) 

I F ( ADEV.GE. 35) CALL 

CKTBKR( 2 , IFIX( ADEV) , CVAL, OF SET, TVAL, 

1 A I M A G ( A I D E V ( IEL, 5) ) ) 
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TVAL = TVAL+2 . 0 

I F ( TVAL2 . LT. TVAL. AND. ND. NE. 0)G0 TO 5 ‘J 1 

TVAL2=TVAL 

ND = 1 

CALL 

DRAWA( 9 00- ( ( TVAL) *180) , (CVAL+FAC1)* 19 0+175.0) 

588 CONTINUE 

******** * **************** * * * * * * * ******** * * * ******* * 

C NOW TO DRAW THE LINE SIDE 

************************************************* 

591 CALL MOVEA(0. 0,0.0) 

AID=(REAL(AIDEV( IEL, 2) ) ) 

IF( AID. GT. 0)CO TO 595 
IDEV=IFIX(-AID) 

STL=REAL(AIDEV(IEL,3)> 

STP=FALT/ 1 0**STL 

IF( REAL(AII)EV(IEL, 5) ) . NE . 0 ) STP=REAL( A I DEV ( IEL , 5) ) 
YM=( STL+FAC)* 190+1 75.0 
XM=0 . 0 

CALL MOVE A( XM , YM ) 

STP=STP* 1 0 
ITP=IFIX( STP) 

D I F= I TP- 1 5 
ID=IFIX(DIF/200.0) 

I F ( ID. LE. 0 ) I D = 1 
TVAL3=10.0 

TD S ET= RE AL ( A I DE V ( IEL, 6) ) 

DO 5 90 K= 1 5 , ITP , II) 

RK = K* . 1 

CALL RELAY ( IDEV , R K , TVAL,TDSET, DEV) 

IF( TVAL. GT. TVAL3 )TVAL=TVAL3 

TVAL3=TVAL 

CALL 

DRAWA( 900- ( ( TVAL+2 .0)*180),((ALOC10(RK)+STL 
1+FAC)* 190+175.0)) 

590 CONTINUE 

IF( REAL ( AIDEV( IEL, 5)).NE.0)CALL 
DRAW A( (900.0) , ( ( A LOG 1 0 ( R K ) + 

1 STL+FAC) *190+175)) 

595 IF(REAL(AIDEV(IEL,2)) .LT.O)GO TO 598 

OF SET = CUR 

ADEV=REAL(AII)EV( IEL, 2 ) ) 

ND = 0 

M=IFIX(80*0FSET) 

KIM=IFIX( 100*AL0G10(FALT) ) 

DO 597 K=M , KIM 

C V A L = K * .01 

I F ( A DEV. LE. 34 ) CALL F U S E ( 1 , AD E V , C V A L , OF SET, TVAL) 
IF(ADEV.GE. 3 5) CAL L 
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CKTBKR( 1 , IFIX(ADEV) , C V AL , OF S ET , T V A L , 

1REAL(AIDEV( I EL , 5) ) ) 

TVAL=TVAL+2 . 0 

IF(TVAL2.LT.TVAL.AND.ND.NE.O)GO TO 601 

TVAL2=TVAL 

ND= 1 

CALL 

DRAWA(900-(TVAL*180) , ( CVAL+FAC) * 1 90+ 1 7 5 . 0 ) 

597 CONTINUE 

601 0FSET=CUR2 

CALL HOVEA(0. 0,0.0) 

NL) = 0 

M=IFIX(80*OFSET) 

DO 599 K=M , KIM 

CVAL=K* .01 

IF( ADEV . LE. 3 4 )CALL FUSE( 2 , ADEV, CVAL, OF SET, TVAL ) 
IF(ADEV.GE.35)CALL 

CKTBKR( 2 , IFIX( ADEV) , CVAL, OF SET, TVAL, 

1REAL( AIDEV ( I EL , 5) ) ) 

TVAL=TVAL+2 . 0 

IF(TVAL2.LT. TVAL. AND. ND.NE.O)CO TO 598 
TVAL2=TVAL 
N D= 1 
CALL 

DRAW A( 9 0 0- ( TVAL* 1 80) , ( CVAL+FAC) * 1 9 0+1 75.0) 

599 CONTINUE 

598 TYPE 600 

600 F0RMAT( ' FINISHED??') 

ND = 0 

CALL MOVEA(0. 0,0.0) 

CALL FLUSH 
ACCEPT 610,1 
610 FORHAT(I) 

IF(I.EQ.O) GO TO 2000 
GO TO 445 

2000 CLOSE( UNIT=30) 

CLOSE ( UNIT=3 1 ) 

CLOSE( U N I T = 3 2 ) 

CLOSE(UNIT=33) 

CLOSE( U N I T = 3 4 ) 

CALL FIHITT(0 , 0) 

STOP 

END 

C SUBROUTINE REFALT ; READS FCUR.DAT TO OUTPUT FAULT 

CURRENTS 

C AT THE BUSSES 

£ * "k & * it * * * * -k * & £ •k it it it * * & it it it it it >’c * it it it it it it it it it it it * it it it it it & it it it 

SUBROUTINE REF ALT( I F ALT , HEL , F ALT , CUR E L ) 
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10 READ( 33 , 1 5) IEF , NN, TFAL, SFAL 

15 FORMAT( 21, 2F) 

1F( IRF.EQ.0)GO TO 20 
IF( IEF.NE. IFALT)G0 TO 10 
I F ( N N . N E . N E L ) G 0 TO 10 
F ALT 1 = TF AL 
FALT2 = S F AL 

F ALT = AMI N 1 ( F ALT 1 , FALT2) 

FL0CUR = 3 . 0* CUREL 

I F ( FALT. LE. FLOCUR) F ALT = AMAX 1 ( F ALT 1 , F ALT2 ) 
FALT=AMAX1 (FALT1 , FALT2) 

20 REWIND 33 

RETURN 
END 



* * k k * * * k k * * k k k k k k k 



k k k kk k k k k k * * * k k k k * k k * * k k k k k k k k k k 



C SUBROUTINE FU S E ; C OM P UTE S TIME VALUE FOR A GIVEN LOAD 

CURRENT 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

SUBROUTINE FU S E ( I D , ADE V , CVAL , OF SET , TVAL) 

TVAL=0 . 0 
IDEV=IFIX( ADEV) 

IF ( IDEV . GT . 1 5 ) GO TO 50 
X= CVAL- 0 F SET 
I F ( ID. NE . 1 )GO TO 40 

TVAL=4 . 0480+(-36 . 4234*X)+( 223 . 7 52 6* (X** 2 ) )+(-75 1 • 279* 

( X**3 ) ) 

l+(1347.01*(X**4))+(-1222.67*(X**5))+(440.813*(X**6>) 
TVAL=TVAL-2 . 0 
GO TO 100 

40 X= C V AL- OF S ET+ . 1 

TVAL=6 . 426+(-37 . 922*X)+( 1 2 1 .909*(X**2 ) )+(-2 1 3 . 420*( X* 

* 3 ) ) + 

1(201 .771* (X**4) >+(-96.983 6*(X**5))+( 18. 6046* (X**6)) 
TV AL= TVA L- 2 . 0 
IF(X.GT. 1 . 4)TVAL=-1 . 95 
GO TO 100 

50 X=CVAL— OFSET 

IF( ID. NE. 1 )GO TO 70 

TVA L= 3 .876 9 7 9 + (-4.262768*X)+(1.21940*(X**2>) 

TVAL= TVA L- 2 . 0 
GO TO 100 

70 X=CVAL-OFSET+. 1 

T V A L = 4 . 614841 + (- 6. 173224 *X)+(3. 2318 06*(X**2)) + (- . 6057 
1 36*(X**3) ) 

TVAL=TVAL-2 . 0 
100 RETURN 

END 

k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 



C SUBROUTINE RELAY ; COMPUTES TIME VALUE FOR A GIVEN 
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CURRENT 

C VALUE FOR A SPECIFIC RELAY 

-k A A * A A A AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



SUBROUTINE RELAY ( NDEV , C U R R , TV A L , T DS ET , DE V ) 
DIMENSION DEV (20,0/4) 

TVAL=0 . 0 
X= ALOG 10(CURR) 

IF(X.GT. 1 . 2)X=1 . 2 
DO 10 1=0,4 

TVAL = TVAL+DEV ( NDEV , I )*(X**I ) 

10 CONTINUE 

IF(NDEV.EQ. 7 )TVAL=TVAL-2 .0 

TVAL=TVAL+TDSET 

RETURN 

END 



AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



C TITLE : LOG. LOG, MAKES LOG LOG FORM FOR GRAPH OF RELAY 

CURVES 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



C OUTLINE BORDER 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



C MUST USE ALOG. LOG, @SYS:GRALIB 



A A A A V? A A AAA A * A A A A A A A * 












SUBROUTINE LOG(IAOFF) 

CALL ERASE 

CALL MOVEA(0. 0,0.0) 

CALL DRAWA(0 .0, 750.0) 

CALL DRAWA( 1000.0, 750.0) 
CALL DRAWA( 1000.0 , 0 . 0) 
CALL DRAWA(0 .0,0.0) 

CALL MOVEA(900. 0,5.0) 

CALL DRAWA( 900. 0,750.0) 

DO 200 N = 0 , 4 
DO 100 M=1 , 10 
DN=1 0* * N 
XM=M*DN 

IF(XM.GT. 35000) GO TO 225 
XX=( ALOGIO(XM) ) * 1 8 0 
XTT=900-XX 

CALL MOVE A( XTT, 175.0) 

CALL DR AW R( 0.0, 575.0) 

100 CONTINUE 

200 CONTINUE 

225 CALL MOVEA(900.0, 175.0) 

CALL DR AW A ( 0.0, 175.0) 

DO 400 11 = 0,2 
DO 30 0 M= 1 , 10 
DN=10**N 
X M = M * D N 
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XX=( ALOCIO(XM) ) * 1 9 O 
XT=XX+1 75 

CALL MOVE A( 95.0, XT ) 

CALL DRAWR(805.0, 0.0) 

300 CONTINUE 

400 CONTINUE 

DO 460 M=0 , 3 
XN=1 75+1 90*M 
CALL MOVEA(900.0,XN) 

CALL NUMBER ( 10, '12' ) 

CALL MOVER(0.0, 10.0) 

J=I AOFF+M 

CALL NUMBER( J , 'll') 

460 CONTINUE 

CALL CHRROT (90) 

CALL MOVEA(960. 0,400.0) 

CALL ANCHOS( ' C1JRRE' , 5 ) 

CALL ANCHOS ( 'NT I ' , 4 ) 

CALL A N C 1 1 0 S ( ' N AMP ',5) 

CALL ANCHOS ( ' S ' , 1 ) 

CALL MOVE A( 9 00. 0,0.0) 

CALL DR AW A ( 900.0,750.0) 

DO 500 M=0 , 5 
XN = 900- 1 3 0 * M 
CALL MOVE A( XN , 150.0) 
IF(M.EQ.O)CALL NUMBER( . 01 , ' F3 . 2 ' ) 
I F ( M . E Q . 1 ) CALL N UM BE R ( . 1 , ' F 2 . 1 ' ) 
IF(M.EQ. 2)CALL NUMBER( 1 , ' 1 1 ' ) 

I F ( M . EQ . 3) CALL N UMBER ( 1 0, ' 12 ' ) 

I F( M . EQ . 4 ) CALL NUMBER ( 100, ' 13' ) 
IF( M. EQ . 5 )C0 TO 510 
500 CONTINUE 

510 CALL MOVEA(20.0, 140.0) 

CALL NUMBER( 1000, ' 14' ) 

CALL CHRROT(O) 

CALL M0VEA(200. 0,755.0) 

CALL ANCHOS ( 'TIME ' , 5) 

CALL ANCHOS ( ' IN SE' , 5) 

CALL ANCHOS( ' CONDS ' , 5 ) 

CALL M0VEA(0. 0,0.0) 

RETURN 

END 






k k k k k k k k 



C SUBROUTINE CKTBKR ; COMPUTES TIME VALUE 

CURRENT 

C VALUE 



k k 

FOR A GIVEN 



kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 



SUBROUTINE CRT BKR ( I D , MD E V , D V AL , OF S ET , T V A L , AM V A L ) 
CVAL=DVAL-OFSET 
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IF(.II). NE. 1 )GO TO 50 
CVAL=CVAL+ .1122 

IF( CVAL.GT. AMVAL. ANI). AllVAL . NE . 0 ) GO TO 90 
I F ( CVAL . LT. 0 ) GO TO 80 

TVA L= 5 . 7 6 90 8 + (-1 0 .882 61 *( CVAL) ) + ( 2 3 . 6 54 1 9 * ( CVA L* * 2 ) ) + 
(-34.7922 

1 * ( CVAL* * 3 ) ) + (33.77472*(CVAL**4)) + (- 21. 1 4 5 3 6 *( CVAL** 5 

) ) + ( 6 . 4 

231778*(CVAL**6)) 

TVAL=TVAL-2 . 0 
GO TO 100 

50 C V A L= C V AL+ . 232 

IF(CVAL.LT.O)GO TO 80 
AMVAL= AMVAL+ .075 

IF( CVAL.GT. AMVAL. AND. AMVAL. NE.O)GO TO 95 

TVAL=5 .85 9 672 2+ (-4. 39 4 3 3 5* CVAL)+ ( - . 50690284* ( CVAL* *2 ) 

)+( 10.40 

1 2 7 3 9 * ( CVAL** 3 ) ) + (- 1 7 . 0 64956* (CVAL* *4 ) ) + ( 1 3. 01 1 54 3* ( C 

VAL**5 ) ) 

2+(-4.010724*(CVAL**6) ) 

TVAL=TVAL-2 . 0 
GO TO 100 
80 TVAL = 6 . 0 

GO TO 150 
90 TVA L=-2 . 0 

GO TO 100 

95 TVAL=ALOC10( .025) 

100 CONTINUE 

150 RETURN 

END 

* * * >V * * * * * & * * & * * -k * * 'k * V* rt * * * * & * * * * * & * * * * * -k * * * * * 



APPENDIX C 



DISK-FILE FORMATS 
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Disk file RES. DAT 



Double precision wire resistance coefficients 



Double precision wire resistance coefficients 



a Double precision wire reactance coefficients 

o 



Double precision wire reactance coefficients 



Disk file TEST . DAT 

element number/line bus/load bus/positive seq p.u. impedance/zero 

seq p.u. impedance/transformer neutral imp line bus/transformer neutral 

imp load bus 

00000 

Load data 

bus number/p. u. VA/per unit subtransient reactance 
0000 

Disk file XFER.DAT 

line bus/load bus/inrush/low withstand/high withstand/ transformer size 
LCUR.DAT 

line bus/load bus/current at line bus/current at load bus 
FCUR.DAT 

fault bus/element number/3 phase fault/1 phase fault 
DEV. DAT 

device number/a /a, /a./a^/a, 
o 1 2 3 A 
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ADEV.DAT 

element number/element bus numbers/device numbers/line bus element 
low pickup, high pickup/load bus Element low pickup, high pickup/ 
magnetic settings line, load/ time dial setting line, load/CT ratios 
line , load/ 



APPENDIX D 



SAMPLE PROGRAM DATA INPUT 
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Before beginning the program sequence, RES. DAT and DEV. DAT 
contain the proper coefficients listed in Table 8.5. All source 
programs have been loaded and saved. PL0TD.F4 uses subroutines con- 
tained in the PDP10 system files so it must be loaded using: 

LOAD PL0TD.F4, @ SYS :GRALIB . 



The first program is started with run INPUT. Underlined words and 
numbers are data input by the user. Words not underlined are prompts 
from the computer in the input sequence that follows. 



Do you want to read from a KYBRD? 

<CR> for input from TTY, 1 for read. . . 
<CR> 



Enter the problem base in MVA 

1 

Enter the problem base in KV 
15 

Enter element data as follows: Bus ( ) to ( ) 

(Length) (Size) (Transformer impedance (P . U. )) (Transformer size 
in MVA) (Transformer High Voltage in KV) (Transformer low 
voltage in KV) 

1 2 3000 6 -(element with no transformer) 

2 3 356 .001 .005 1 15 4.1 -(element with a transformer) 

Input transformer ZN in P.U. refer to BUS 

Side, (BUS 2) (BUS 3) . . . 

.001 .005 .001 .005 



(This is continued for each element until all elements of 
the type in RES. DAT have been entered) 

0^ -(tells computer no more elements) 

At this point if another wire type is to be used the 
following procedure applies if not skip this section 
and proceed to (l) 
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Input Bus 1 voltage. . .0^ 

Data is accepted and a list follows 

(Here the computer types all computed information on the 
elements entered so that they may be checked for error.) 

Now enter the load data in the following manner;. . . 

(Bus) (Load in KVA) (Pwr Factor). . .or. . . (Bus) (Load in HP) 

(Pwr Factor). . .or. . . (Bus) (Load in Amps) (Pwr Factor) 

What method will you use?. . .VA, HP, or Amps? 

S_ -(writes element data computed to disk) 

End of execution 

(Load RES. DAT with the new wire type coefficients) 

RUN INPUT 

Do you want to read from a KYBRD? 

<CR> for input from a TTY, 1 for read. . . 

1 _ 

Old data is entered ready for new line data 
Enter element data as follows. . .(Program runs as 
in the first section) 

Picking up the program from completion of element data, 
the data entry follows. 

Input Bus 1 Voltage. . . 15QQQ 
Input Bus 2 Voltage. . . 15Q0Q 
Input Bus 3 Voltage. . . 4100 

Input 'Bus 8 Voltage. . . 440 

Data is accepted and a list follows 

(Computer types computed element data) . 

Number of buses is 8 

Now enter the load data in the following manner;. . . 

(Bus) (Load in KVA) (Pwr Factor). . .or. . . (Bus) (Load in HP) 

(Pwr Factor). . .or. . . (Bus) (Load in Amps) (Pwr factor) 

What method will you use? KVA, HP, or Amps? 

HP 

Enter. . . 

4 800 .80 

What method will you use?. . .KVA, HP, or Amps? 

HP 

8 500 .90 

(This is continued until all loads are entered? 

What method will you use?. . .KVA, HP, or Amps? 

End of Execution 
RUN LDFLO 
End of Execution 
RUN ZBUS 

Enter the zero sequence impedance for the loads at each bus. . . 
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8 .001 .005 -for each load where the zero sequence impedance 

• will not be assumed to equal the positive sequence 

impedance 

0 -tells the computer there are no more entries 

Enter the mutually coupled elements in per unit in the following 
format, TEST.DAT has the element numbers. . (Element) (Element ) 
(Zero Sequence Value) 

230. 001 (The zero must be entered because the coupling is a 
reactive value) 

0 _ -(no more or no coupling) 

End of execution 
RUN FALT 

Enter the fault impedance (Real) (Imaginary) 

0 0 

End of execution 
RUN COORD 

All load currents are ready, use device number list 
Select device for the load at bus 4 
4 



Select the device for the load at Bus 8 

2 

Input temperature for ambient ahift, 1 for less than 
110 Deg C, 2 for 60 Deg C 
<CR> for no shift. . . 

<CR> 

Does the device have a magnetic element ? (0) Yes 
0 

Select device for Bus 6 protecting element feeding load at Bus 8 

2 

Select the Westinghouse Co. relay, refer to Table 2 
Does the relay have an inst. element if yes 0 
0 



(Until all backup devices have been set) . 
Enter the device for the Bus 6 of element 5 
_4 

Select the Westinghouse Co. relay, refer to Table 2 
Does the relay have an inst. element if yes 0 
0 



(Until all devices in the system have been set). 

Device listing -(Here the computer lists the matrix ADEV) 

End of execution 
Run PLOTD 

Devices have been set if you want to plot graph paper 0 if not 
type 1 
0 
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This will plot curves based on the lowest voltage in the 
system 

Input the current offset for the graph, integer 1_ 

Input element number 

1_ 

Finished? 

1_ -(no, 0 yes) 

Devices have been set if you want to plot graph paper 

0 if not type 1_ 

1 

Input element number 

6 

Finished? 

0 

End of execution. 

If a device did not plot as desired its settings can be changed 
in ADEV.DAT. 
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